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Preface 


Differential equations belong to one of the 
main mathematical concepts. They are 
equations for finding functions whose deriv- 
atives (or differentials) satisfy given condi- 
tions. The differential equations arrived at 
in the process of studying a real phenomenon 
or process are called the differential model 
of this phenomenon or process. It is clear 
that differential models constitute a particu- 
lar case of the numerous mathematical models 
that can be built as a result of studies of 
the world that surrounds us. It must be 
emphasized that there are different types of 
differential models. This book considers none 
but models described by what is known 
as. ordinaru, differential eayations. one. char- 
acteristic of which is that the unknown 
functions in these equations depend on 
a single variable. 

In constructing ordinary differential mo- 
dels it is important to know the laws of 
the branch of science relating to the nature 
of the problem being studied. For instance, 
in mechanics these may be Newton’s laws, 
in the theory of electric circuits Kirchhoff's 
laws, in the theory of chemical reaction 
rates the law of mass action. 

Of course, in practical life we offen have 
to deal with cases where the laws that en- 
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able building a differentia lequation (or sev- 
eral differential equations) are not known, 
and we must resort to various assumptions 
(hypotheses) concerning the course of the 
process at small variations of the parame- 
ters, the variables. Passage to the limit 
will then lead to a differential equation, 
and if it so happens that the results of 
investigation of the differential equation as 
the mathematical model agree with the ex- 
perimental data, this will mean that the 
hypothesis underlying the model reflects 
the true situation.* 

When working on this book, I had two 
goals in mind. The first was to use examples, 
rich in content rather than purely illustra- 
tive, taken from various fields of knowledge 
so as to demonstrate the possibilities of 
using ordinary differential equations in 
gaining an understanding of the world about 
us. Of course, the examples far from ex- 
haust the scope of problems solvable by or- 
dinary differential equations. They give 
an idea of the role that ordinary differential 
equations play in solving practical problems. 

The second goal was to acquaint the read- 


* If the reader wishes to know more about mathe- 
matical models, he can turn to the fascinating books 
by A.N. Tikhonov and D.P. Kostomarov, Stories 
About Applied Mathematics (Moscow: Nauka, 
1979) (in Russian), and N.N. Moiseev, Mathemat- 
ics Stages an Experiment (Moscow: Nauka, 1979) 
(in Russian), 
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er with the simplest tools and methods 
used in studying ordinary differential equa- 
tions and characteristic of the qualitative 
theory of differential equations. The fact is 
that only in rare cases are we able to solve 
a differential equation in the so-called 
closed form, that is, represent the solution 
as a formula that employs a finite number 
of the simplest operations involving ele- 
mentary functions, even when it is known 
that the differential equation has a solu- 
tion. In other words, we can say that the 
great variety of solutions to differential 
equations is such that for their representa- 
tion in closed form a finite number of ana- 
lytical operations is insufficient. A simi- 
lar situation exists in the theory of algebraic 
equations: while for first- and second-order 
algebraic equations the solutions can always 
be easily expressed in terms of radicals 
and for third- and fourth-order equations 
the solutions can still be expressed in terms 
of radicals (although the formulas become 
very complicated), for a general algebraic 
equation of an order higher than the fourth 
the solution cannot be expressed in radicals. 

To return to differential equations. If 
an infinite series of this or that form is used 
to represent the solutions, then the scope 
of solvable equations broadens considerab- 
ly. Unfortunately, it often happens that 
the most essential and interesting proper- 
ties of the solutions cannot be revealed by 
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studying the form of the series. More than 
that, even if a differential equation can be 
solved in closed form, more often than not 
it is impossible to analyze such a solution 
since the relationship between the various 
parameters of the solution often proves 
to be extremely complicated. 

This shows how important it is to develop 
methods that make it possible to acquire 
the data on the various properties of the 
solutions without solving the differential 
equations themselves. And indeed, such 
methods do exist. They constitute the es- 
sence of the qualitative theory of differen- 
tial equations, which is based on the gener- 
al theorems regarding the existence and 
uniqueness of solutions and the continuous 
dependence of solutions on the initial data 
and parameters. The role of existence and 
uniqueness theorems is partially discussed 
in Section 2.2. As for the general qualitative 
theory of ordinary differential equations, 
ever since J.H. Poincaré and A.M. Lya- 
punov laid the foundations at the end of 
the 19th century, the theory has been inten- 
sively developing and its methods are wide- 
ly used when studying the world about us. 

I am indebted to Professors Yu.S. Bog- 
danov and M.V. Fedoryuk for the construc- 
tive remarks and comments expressed in 
the process of preparing the book for publi- 
cation, 

VV. Amel’kin 


Chapter 1 


Construction 
of Differential Models 
and Their Solutions 


1.1 Whose Coffee Was Hotter? 


When Tom and Dick ordered coffee and 
cream in a lunch room, they were given both 
simultaneously and proceeded as follows. 
Tom poured some of his cream into the cof- 
fee, covered the cup with a paper napkin, 
and went to make a phone call. Dick cov- 
ered his cup witha napkin and poured the 
same amount of cream into his coffee only 
after 10 minutes, when Tom returned. The 
two started drinking their coffee at the same 
time. Whose was hotter? 

We willsolve this problem on the natural 
assumption that according to the laws of 
physics heat transfer through the surface of 
the table and the paper napkin is much less 
than through the sides of the cups and that 
the temperature of the vapor above the sur- 
face of the coffee in the cups equals the tem- 
perature of the coffee. 

We start by deriving a relationship indi- 
cating the time dependence of the tempera- 
ture of the coffee in Dick’s cup before the 
cream is added. 

In accordance with our assumption on 
the basis of a law of physics, the amount 
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of heat transferred to the air from Dick’s 
cup is determined by the formula 

dQ = nt sat, (1) 
where 7 is the coffee’s temperature at time 
t, 6 the temperature of the air in the lunch 
room, y the thermal conductivity of the 
material of the cup, J the thickness of the 
cup, and s the area of the cup’s lateral sur- 
face. The amount of heat given off by the 
coffee is 


dQ = —cm dT, (2) 


where c is the specific heat capacity of the 
coffee, and m the mass (or amount) of coffee 
in the cup. If we now consider Eqs. (1) 
and (2) together, we arrive at the following 
equation: 


n a sdt= —cmdT, 


which, after variable separation, can be re- 
written as follows: 
dT. ns 

T—6 — Lem df. (3) 
Denoting the initial temperature of the coffee 
by 7, and integrating the differential equa- 
tion (3), we find that 
T=0+(T,—6) exp (——* t). (4) 


lem 


This formula is the analytical expression 
of the law whereby the temperature of the 
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coffee in Dick’s cup varied prior to the addi- 
tion of the cream. 

WB Now let us establish the law whereby the 
temperature of the coffee in Dick’s cup 
changed after Dick poured in cream. We 
use the heat balance equation, which here 
can be written as 


em (T — 8p) = cm, (8p — T)), (9) 


where 0p is the temperature of the Dick’s 
coffee with cream at time ¢, 7’, the temper- 
ature of the cream, c, the specific heat ca- 
pacity of the cream, and m, the mass of 
cream added to the coffee. 

Equation (5) yields 


_— yy 
(i= taLam tr aoe (6) 


Bearing in mind formula (4), we can re- 
write (6) as follows: 


Cymy, 
Op =P, 


em cm, Ee Tl 
x | 8+ (T>— 8) exp( — t)], (7) 


which constitutes the law whereby the tem- 
perature of the coffee in Dick’s cup varies 
after cream is added. 

To derive the law for temperature varia- 
tion of the coffee in Tom’s cup we again 
employ the heat balance equation, which 
now assumes the form 


cm (Ty — 89) = cym, (8) — 71), (8) 
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where 0, is the temperature of the mixture. 
If we solve (8) for 0), we get 


cm, cm 

80 om-Feym, 11 em Fem, 0” 

Then, using Eq. (4) with 0, serving as the 
initial temperature and cm + c,m, substi- 
tuted for em, we arrive at the law for the 
temperature (@;) variation of the coffee 
in Tom’s cup as analytically given by the 
following formula: 


Op = 0+ | 7, +" _ 7,8 


cm-+cym, cm—+cym, 
Segre | ieee 
x el ( l (cm + C4mM)) t} : (9) 


Thus, to answer the question posed in 
the problem we need only turn to formulas 
(7) and (9) and carry out the necessary cal- 
culations, bearing in mind that c, ~ 3.9 x 
10? J/kg-K, c~ 4.1 x 10® J/kg-K, and 
n 0.6 V/m-K and assuming, for the sake 
of definiteness, that m, = 2 x 10~* kg, 
m=8 X 10-* kg, T, = 20°C, 0 = 20 C, 
Ty = 00°C, sit x 107° m*,. and 7 = 
2x 10-3 m. The calculations show that 
Tom’s coffee was hotter. 


1.2 Steady-State Heat Flow 


The reader will recall that a steady-state 
heat flow is one in which the object’s tem- 


= 
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Fig. 1 Fig. 2 


perature at each point does not vary with 
time. 

In problems whose physical content is 
related to the effects of heat flows, an impor- 
tant role is played by the so-called isother- 
mal surfaces. To clarify this statement, let 
us consider a heat-conducting pipe (Figure 
1) 20 cm in diameter, made of a homogene- 
ous material, and protected by a layer of 
magnesium oxide 10 cm thick. We assume 
that the temperature of the pipe is 160 °C 
and the outer surface of the protective cov- 
ering has a temperature of 30°C. It is 
intuitively clear then that there is a surface, 
designated by a dashed curve in Figure 
2, at each point of which the temperature 
is the same, say 95 °C. The dashed curve in 
Figure 2 is known as an isotherm, while the 
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surface corresponding to this curve,is known 
as an isothermal surface. In general, iso- 
therms may have various shapes, depend- 
ing, for one, on the nonsteady-state nature 
of the heatfflow and on the nonhomogeneity 
of the material. In the case at hand the 
isotherms (isothermal surfaces) are represent- 
ed by concentric circles (cylinders). 

We wish to derive the law of temperature 
distribution inside the protective coating 
and find the amount of heat released by the 
pipe over a section 1 m long in the course 
of 24 hours, assuming that the thermal con- 
ductivity coefficient k is equal to 1.7 x 
10-4, 

To this end we turn to the Fourier law, 
according to which the amount of heat re- 
leased per unit time by an object that is in sta- 
ble thermal state and whose temperature T at 
each point is solely a function of coordinate x 
can be found according to the formula 


= —kF (z) -. = const, (10) 


where F (x) is the cross-sectional area normal 
to the direction of heat flow, and k the thermal 
conductivity coefficient. 

The statement of the problem implies 
that F (x) = 2nzl, where | is the length of 
the pipe (cm), and z the radius of the base 
of the cylindrical surface lying inside the 
outer cylinder. Then on the basis of (10) 
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we get 
30 0 20 
oo 
\ df = ~ 0.00017 x 2a a (11) 
160 10 
T 0 x ae 
| dt= —oypirxuar \ E- 7) 
160 10 


Integrating (11) and (12) yields 


1460—T InO1x _ log0.iz 
130 ~ In2 ~~ Jog2 


Hence 
T = 591.8 — 431.8 log z. 


This formula expresses the law of tempera- 
ture distribution inside the protective 
coating. We see that the length of the pipe 
plays no role in this law. 

To answer the second question, we turn 
to Eq. (11). Then for / = 100 cm we find 
that 


__ 130 x 0.00017 X 2 x 100 


Q In 2 


200 x 130 x 0.00017 
=" ~"0.69315 


and, hence, the amount of heat released by 
the pipe in the course of 24 hours is 24 X 
60 x 602 = 726 852 J. 


2—0770 
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1.3 An Incident in a National Park 


While patrolling a national park, two for- 
est rangers found the carcass of a slain wild 
boar. Inspection showed that the poacher’s 
shot had been precise and the boar had died 
on the spot. Reasoning that the poacher 
should return for the kill, the rangers de- 
cided to wait for him and hid nearby. Soon 
two men appeared, heading directly for the 
dead boar. When confronted by the rangers, 
the men denied having anything to do with 
the poaching. But by this time the rangers 
had collected indirect evidence of their 
guilt. It only remained to establish the 
exact time of the kill. This they did by 
using the law of heat emission. Let us see 
what reasoning this involved. 

According to the law of heat emission, 
the rate at which an object cools off in air 
is proportional to the difference between 
the temperature of the object and that of 
the air, 


= —k(x—a), (13) 
where z is the temperature of the object at 
time ¢, a the temperature of the air, and 
k a positive proportionality coefficient. 

Solution of the problem lies in an analy- 
sis of the relationship that results from in- 
tegrating the differential equation (13). 
Here one must bear in mind that after the 
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boar was killed, the temperature of the air 
could have remained constant but also 
could have varied with time. In the first case 
integration of the differential equation (13) 
with variables separable leads to 


t—a 


In 


= —kt, ra, (14) 


where x, is the temperature of the object 
at time t = 0. Then, if at the time when 
the strangers were confronted by the rangers 
the temperature of the carcass, x, was 34 °C 
and after an hour 29 °C and if when the shot 
was fired the temperature of the boar was 
x = 37°C and the temperature of the air 
a = 214 °C, we can establish when the shot 
was fired by putting t = 0 as the time when 
the strangers were detained. Using these 
data and Eq. (14), we find that 


31 — 21 
29— 24 


A= In = In 1.25 = 0.22314. (15) 


Now, substituting this value of k and z = 
o7 into (44), we get 


t= — 59314 ‘8 3p aT 
‘ 
= ~ 959314, In 1.6= — 2.10630. 


In other words, roughly 2 hours and 6 min- 
utes passed between the time the boar was 
killed and the time the strangers were de- 
tained. 


2" 
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In the second case, that is, when the tem- 
perature of the air varies in time, the cool- 
ing off of the carcass is expressed by the 
following nonhomogeneous linear differen- 
tial equation 


S 4 ke = ka(t), (16) 


where a (t) is the temperature of the air at 
time 7. 

To illustrate one of the methods for deter- 
mining the time when the boar was killed, 
let us assume that the temperature of the 
carcass was 30 °C when the strangers were 
detained. Let us also suppose that it is 
known that on the day of the kill the tem- 
perature of the air dropped by 4 °C every 
hour in the afternoon and was 0 °C when 
the carcass was discovered. We will also 
assume that after an hour had passed after 
the discovery the temperature of the carcass 
was 20°C and that of the air was down 
to —1° C. If we now assume that the shot 
was fired at ¢=0O and that x) = 37°C 
at ¢=0, we get a(t) = t* —t, where 
t = t* is the time when the rangers dis- 
covered the carcass. 

Integrating Eq. (16), we get 


ga (87 —~t—kYe* 4 e* —¢ + kL 


If we now bear in mind that x = 30°C 
at ¢=¢* and «x = 25°C at t= t* +1, 
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the last formula yields 


(37 — t* — k-) exp (—kt*) +k? = 30, 
(37 — t* — k-') exp [—k& (t* + 1)] 
+ i-l = 26. 


These two equations can be used to derive 
an equation for k, namely, 


(30 — k-e-* — 26 +k = 0. (17) 


We can arrive at the same equation start- 
ing from different assumptions. Indeed, 
suppose that at ¢=0 the carcass of the slain 
boar was found. Then a(t) = —t and 
we arrive at the differential equation 


St 4 k= — kt (18) 


(with the initial data z), = 30 at t = 0), 
from which we must find x as an explicit 
function of ¢. 

Solving Eq. (18), we get 


x = (30 — k-') exp (—At) —t + k-!. (49) 


Setting « = 25 and ¢ = 1 in this relation- 
ship, we again arrive at Eq. (17), which ena- 
bles solving the initial problem numerical- 
ly. 
Indeed, as is known, Eq. (17) cannot be 
solved algebraically for k. But it is easily 
solved by numerical methods for finding the 
roots of transcendental equations, for one, 
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by Newton’s method of approximation. This 
method, as well as other methods of succes- 
sive approximations, is a way of using a 
rough estimate of the true value of a root 
used to obtain more exact values of the 
root. The process can be continued until 
the required accuracy is achieved. 

To show how Newton’s method is used, 
we transform Eq. (17) to the form 


30k — 41 + (1 — 26k) exp(k) =0, — (20) 


and Eq. (19), setting x = 37, to the form 

(87k — 1 + kt) exp (kt) — 30k +1 = 0. 
(21) 

Both equations, (20) and (21), are of the 

type 

(ax + b) exp (Av) + cr +d = 0. (22 


) 
If we denote the left-hand side of Eq. (22) 
by @ (z), differentiation with respect to z 
yields 
p (x) = (Aax + db + a) exp (Ax) + ¢, 
p” (xz) = (ax + 2b + 2da) exp (Az). 
Then according to Newton’s method for 


finding a root of Eq. (22), if for the ith 
approximation z; we have the inequality 


P (zi) P" (x) > 0, 
the next approximation, x;+,, can be found 
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via the formula 


__ 9 (zi) 

p’ (xj) ° 

To directly calculate the root (to within, 
say, one part in a million), we compile 
the following program using BASIC*: 


10 CLS:PRINT “Solution of equation by New- 
ton’s method” 

45 INPUT “lambda=”, L; 

20 INPUT “a=”, A: INPUT “b=”, B; 

30 INPUT “c=”, C: INPUT “d=”, D 

40 INPUT “approximate value of root=”, X 

50 PRINT “xX”, “f", “f’", “f’" 

100 E = EXP (L*X) 

110 F = (A*X + B)*¥E+ C*X + D 

1420 Fi = (L*¥(A*X -+ B) + A)*E4+C 

130 F2 = L* (L* (A*X + B) -+ 2*A)*E 

150 PRINT X, F, Fi, F2 

151 IF F = 0 THEN END 

155 IF Fi = 0 THEN PRINT “Newton’s method 
is divergent”: END 

170 IF F*F2 <0 THEN PRINT “Newton's 
method is divergent”: END 

190 X = X — F/F1 

200 GOTO 100 


Lisy = LT; 


In this program X, f, f’ and f” stand for 
kn» © (kn), @ (kp), and p” (k,), respective- 
ly. 

* Editor's note. Some BASIC lines, e.g. line 


10, have been split due to the printed format; 
they should be input as one line. 
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After starting the program enter the re- 
quested values of the coefficients of the equa- 
tion and the initial value of the root. The 
results can be listed in a table of the approx- 
imate values of the root and the respective 
values of the function and its first and 
second derivatives. 

_ Employing this general procedure, let us 
turn to Eq. (20). Differentiating its left-hand 
side @ (k) with respect to k, we arrive at 


yp’ (k) = 30 — (25 + 26k) exp (k). 


It can then easily be verified that «¢ (0) = 
0, ~(1)< 0, and q’ (0) >0. Thus, the 
function @ increases in a small neighborhood 
of the origin and then decreases to a nega- 
tive value at k =1.Thisimplies that in the 
interval (0, 1) there is a root of the equation 
@p (k) = 0. To find this root we run the pro- 
gram. Below we give the protocol of the 
program: 


Solution of equation by Newton’s method 
lambda= 14 


a=—26 b=1c=30d=—1 
Approximate value of root = 0.5 

X f i f” 
5 —5.784655 —32.65141 —105.5182 
.322836 —1.525956 —16.11805 —82.02506 
.2281622 —.3514252 —8.859806 —71.52333 


.1884971 —5.519438E—02 —6.103379 —67.49665 
4794939 —e.748013E—03 —2,497021 —66,60768 
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178954 —7.629395KE—06 —5.46373 —66.55884 
.1789526 —4.768372E—07 —5.463642 —66.5587 
1789525 =O —5.463635 —66.5587 
OK 


In this protocol X, f, f’, and f” stand for 
kn» @ (kn), 9 (kn), and op” (k,), respective- 
] 


The final step in solving the problem con- 
sists in substituting the calculated value 
hk, wk ~ 0.178952 into Eq. (24) and solv- 
ing the latter for ¢ (the time when the wild 
boar was killed). To employ the above 
scheme, we denote the left-hand side of Eq. 
(21) by g(t). Then, selecting —1 for the 
value of ¢, and bearing in mind that in 
this case a = k,b = 37k —1,c =0,d = 
—30k+1,and 4=k, we can find the time 
when the boar was killed. To this end 
we find the coefficients b and d, which prove 
to be equal to 5.621243 and —4.368975, re- 
spectively, and then running the above pro- 
gram, find the sought time t. The proto- 
col of the program is given below: 


Solution of equation by Newton’s method 


lambda = 0.1789525 
a = 0.1789525 b = 5.621243 c=0 


d = —4.368575 
Approximate value of root = —4 

xX f f’ i” 
—1 181972 9639622 .1992802 


1.188775 3.506184E—03 .9270549 .1917861 
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Fig. 3 

—1.192557 1.430512E—06 .9263298  .1916388 
—1.192559 0 9263295 .1916387 
OK 


In this protocol X, f, f’, and f£” stand for 
tr, & (tn), g’ (tn), and g” (t,), respectively. 
These results imply that the boar was 
killed approximately 4 hour and 12 minutes 
before the rangers discovered the carcass. 


1.4 Liquid Flow Out of Vessels. 
The Water Clock 


The two problems that we now discuss illu- 
strate the relationship between the physical 
content of a problem and geometry. But first 
let us examine some general theoretical con- 
clusions. 

We take a vessel (Figure 3) whose hori- 
zontal cross section has an area that is a 
function of the distance from the bottom 
of the vessel to the cross section, Suppose 
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that initially at time ¢ = 0 the level of the 
liquid in the vessel is at a height of h 
meters. We will also suppose that the area 
of the vessel’s cross section at height zx 
is denoted by S (x) and that the area of the 
opening in the bottom of the vessel is s. 

As is known, the rate v at which the liquid 
flows out of the vessel at the moment when the 
liquid’s level is at height x is given by the for- 
mula v=kY2gz, where g =9.8 m/s?, 
and k is the rate constant of the outflow pro- 
cess. 

In the course of an infinitesimal time 
interval dé the outflow of the liquid can be 
assumed uniform, whereby during dt a 
column of liquid with a height of vdt 
and a cross-sectional area of s will flow out 
of the vessel, which causes the level of the 
liquid to change by —dz (the “minus” be- 
cause the level lowers). 

The above reasoning leads us to the fol- 
lowing differential equation 


ks V 2gx dt = —S (z) dz, 
which can be rewritten as 
S (z) 


di SS ee 
ks Y 2gx 


(23) 


Let us now solve the following problem. 
A cylindrical vessel with a vertical axis 
six meters high and four meters in diameter 
has a circular opening in the bottom. The 
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radius of this opening is 1/12 m. Find how 
the level of water in the vessel depends on 
time ¢ and the time it takes all the water to 
flow out. 

By hypothesis, S (7) =4n and s= 
1/144. Since for water k = 0.6, Eq. (23) 
assumes the form 


di = acon dz. 
Vz 


Integrating this differential equation yields 
t=434.304[V6—Va], O0<x<6, 


which is the sought dependence of the level 
of water in the vessel on time ¢. If we put 

= (Q in the last formula, we find that it 
takes approximately 18 minutes for all the 
water to flow out of the vessel. 

Now a second problem. An ancient water 
clock consists of a bowl with a small hole 
in the bottom through which water flows 
out of the bowl (Figure 4). Such clocks were 
used in ancient Greek and Roman courts to 
time the lawyers’ speeches, so as to avoid 
prolonged speeches. We wish to determine 
the shape of the water clock that would en- 
sure that the water level lowered at a con- 
stant rate. 

This problem can easily be solved via 
Eq. (23). We rewrite this equation in the 
form 


C2, dz 


= 
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Fig. 4 


Precisely, if we suppose that the bow! has 
the shape of a surface of revolution, in ac- 
cordance with the notations used in Figure 
4, Eq. (24) yields the following result: 


Vi = oe aie a, (25) 


where a = v, = dz/dt is the projection on 
the x axis of the rate of motion of the 
water's free surface, which is constant by 
hypothesis. Squaring both sides of Kq. 
(25), we arrive at the equation 


2 = cr’, (26) 


with c = a*n?/2gk?s*. The latter means that 
the sought shape of the water clock is ob- 
tained by rotating curve (26) about the zx 
axis. 
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1.5 Effectiveness of Advertising 


Suppose that a retail chain is selling com- 
modities of a certain type, say B, about 
which only z buyers out of N potential buy- 
ers know at time ¢. Let us also assume that 
to speed up sales the chain has placed pro- 
motion materials in the local TV and radio 
network. All further information about the 
commodities is distributed among the buy- 
ers via personal contact between them. We 
may assume with a high probability that 
after the TV and radio network have re- 
leased the information about B, the rate of 
change in the number of persons knowing 
about B is proportional both to the number 
of buyers knowing about B and to the num- 
ber of buyers not knowing about B. 

If we suppose that time is reckoned from 
the moment when the promotion materials 
and advertisements were released and 
N/y persons have learned about the commo- 
dities, we arrive at the following differen- 
tial equation: 


<* —kx (N— 2), (27) 


with the initial condition that x = N/y 
at t = 0. In this equation k is a positive 
proportionality factor. Integrating, we find 
that 

4 x 


W In = =khi+c. 


= 
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Assuming that NC = C,, we arrive at the 
equation 


zx 
— ApNht 
Noe AeNht, 


with A =e. When solved for z, this 
equation yields 


Nht 
aN cee (28) 


with P = 4/A. 

In economics, Eq. (28) is commonly 
known as the logistic-curve equation. 

If we allow for the initial condition, Eq. 
(28) becomes 


Se ae 
1 (y—1) eo NMFS 


Figure 5 provides a schematic of a logis- 
tic curve for y = 2. In conclusion we note 
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that the problem of dissemination of tech- 
nological innovations can also be reduced 
to Eq. (27). 


1.6 Supply and Demand 


As is known, supply and demand constitute 
economic categories in a commodity econo- 
my, categories that emerge and function on 
the market, that is, in the sphere of trade. 
The first category represents the commodi- 
ties that exist on the market or can be deliv- 
ered to it, while the second represents the 
demand for commodities on the market. 
One of the main laws of a market economy 
is the law of supply and demand, which 
can be formulated simply by saying that 
for each commodity some price must exist 
that will cause the supply and the demand 
to be just equal. Such a price establishes an 
“equilibrium” on the market. 

Let us consider the following problem. 
Suppose that in the course of a (relatively 
long) time interval a farmer sells his produce 
(say, apples) on the market, and that he 
does this immediately after the apple har- 
vest has been taken in and then once each 
following week (i.e. with a week’s interval). 
The farmer’s stock of apples being fixed af- 
ter he has collected his harvest, the week’s 
supply will depend on both the expected 
price of apples in the coming week and the 
expected change in price in the following 
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weeks. If it is assumed that in the coming 
week the price of apples will fall and in the 
following weeks it will grow, then supply 
will be restrained provided that the ex- 
pected rise in price exceeds storage costs. In 
these conditions the greater the expected rise 
in price in the following weeks, the lower 
the supply of apples on the market. On the 
other hand, if the price of apples in the 
coming week is high and then a fall is ex- 
pected in the following weeks, then the 
greater the expected fall in price in the fol- 
lowing weeks, the higher the supply of ap- 
ples on the market in the coming week. 
If we denote the price of apples in the 
coming week by p and the time derivative 
of price (the tendency of price formation) 
by p’, both supply and demand are func- 
tions of these quantities. As practice has 
shown, depending on various factors the 
supply and demand of a commodity may be 
represented by different functions of price 
and tendency of price formation. For exam- 
ple, one function is given by a linear de- 
pendence expressed mathematically in the 
form y = ap’ + bp +c, where a, b, and 
ec are real numbers (constants). Say, in 
our example, the price of apples in the 
coming week is taken at 1 rouble for 4 kg, 
in t weeks it was p (t) roubles for 1 kg, 
and the supply s and demand q are given 
by the functions 
s=44p' + 2p—1, q=4p' — 2p4 39. 
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Then, for the supply and demand to be in 
equilibrium, we must require that 


44p' + 2p —1 = 4p’ — 2p 4 39. 


This condition leads us to the differential 
equation 


dp 
Integration yields p = Ce-!% + 10. If we 
allow for the initial condition that p = 1 
at t = 0, the equilibrium price is given by 
the formula 


p = —9e-1! + 10, (29) 


Thus, if we want the supply and demand to 
be in equilibrium all the time, the price 
must vary according to formula (29). 


1.7 Chemical Reactions 


A chemical equation shows how the inter- 
action of substances produces a new sub- 
stance. Take, for instance, the equation 


oH, + O, > 2H,0. 


It shows how as a result of the interaction 
of two molecules of hydrogen and one mo- 
lecule of oxygen two molecules of water are 
formed. 
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Generally a chemical reaction can be writ- 
ten in the form 


aA + 6B + cC -+- +>mM + nN 

+ pP + 

where A,B,C, are molecules of the in- 
teracting substances, M, N, P, .. mole- 
cules of the reaction products, and the con- 
stants a, b, c, .., m,n,p,. positive 


integers that stand for the number of mole- 
cules participating in the reaction. 

The rate at which a new substance is 
formed in a reaction is called the reaction 
rate, and the active mass or concentration 
of a reacting substance is given by the num- 
ber of moles of this substance per unit vol- 
ume. 

One of the basic laws of the theory of 
chemical reaction rates isthe law of mass 
action, according to which the rate of a chem- 
ical reaction proceeding at a constant tem- 
perature is proportional to the product of the 
concentrations of the substances taking part 
in the reaction at a given moment. 

Let us solve the following problem. Two 
liquid chemical substances A and B occu- 
pying a volume of 40 and 20 liters, respec- 
tively, form as a result of a chemical reaction 
anew liquid chemical substance C. Assum- 
ing that the temperature does not change 
in the process of the reaction and that every 
two volumes of substance A and one vol- 
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ume of substance B form three volumes of 
substance C, find the amount of substance 
C at an arbitrary moment of time ¢ if it 
takes 20 min to form 6 liters of C. 

Let us denote by x the volume (in liters) 
of substance C that has formed by the time 
t (in hours). By hypothesis, by that time 
22z/3 liters of substance A and 2/3 liters of 
substance B will have reacted. This means 
that we are left with 10 — 22/3 liters of A 
and 20 — 2/3 liters of B. Thus, in accord- 
ance with the law of mass action, we arrive 
at the following differential equation 


dz 22x x 
$= (10-4) (20-§) 
which can be rewritten as 
d 

Fr = (15 —2) (60—2), 


where & is the proportionality factor (k = 
2K/9). We must also bear in mind that since 
initially (¢ = 0) there was no substance 
C, we may assume that z = 0 att = O. As 
for the moment ¢t = 1/3, we have x = 6. 

Thus, the solution of the initial problem 
has been reduced to the solution of a so- 
called boundary value problem: 


 —k(15—z)(60—z), 2(0)=0, 
x (4/3) = 6. 
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To solve it, we first integrate the differen- 
tial equation and allow for the initial con- 
dition z (0) =O. As a result we arrive 
at the relationship 


60—-z 
=s Abdkt 
oe == 4eAokt, 


Since zc = 6 at t = 1/3, substituting these 
values into the relationship yields e1®* = 
30/2. Hence, 


60 — zx 
fea — 4 (e15ky3t . 4 (3/2)3t, 
that is, 


1 — (2/3)%! 
1 (1/4) (2/3)34 


This gives the amount of substance C 
that has been formed in the reaction by 
time tf. 

A remark is in order. From practical 
considerations it is clear that only a finite 
volume of substance C can be formed when 
10 liters of A interact chemically with 20 
liters of B. However, a formal study of 
the above relationship between zx and ¢ 
shows that for a finite f¢, namely at (2/3)3! = 
4, the variable x becomes infinite. But this 
fact does not contradict practical consid- 
erations because it is realized only for a 
negative value of t, while the chemical reac- 
tion is considered only for ¢ nonnegative. 
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1.8 Differential Models in Ecology 


Ecology studies the interaction of man and, 
in general, living organisms with the envi- 
ronment. The basic object in ecology is the 
evolution of populations. Below we describe 
differential models of populations that deal 
with their reproduction or extinction and 
with the coexistence of various species of 
animals in the predator vs. prey situation.* 

Let z(t) be the number of individuals 
in a population at time ¢. If A is the number 
of individuals in the population that are 
born per unit time and B the number of 
individuals that die off per unit time, then 
there is reason to say that the rate at which 
x varies in time is given by the formula 


d 

= = AB. (30) 
The problem consists in finding the depen- 
dence of A and B on x. The simplest situa- 
tion is the one in which 


A=az, B = bz, (31) 
where a and Dare the coefficients of births 
and deaths of individuals per unit time, 


respectively. Allowing for (31), we can re- 
write the differential equation (80) as 


* For more detail see J.D. Murray, “Some simple 
mathematical models in ecology,” Math. Spectrum 
16, No. 2: 48-54 (1983/84). 
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Assuming that at ¢ = t, the number of in- 
dividuals in the population is z =z, 
and solving Eq. (32), we get 


x (t) = x, exp [(a — b) (t — t,)]. 


This implies that if a > 6, then the number 
of individuals zx tends to infinity as t > oo, 
but if a< b, then z—>0O as t- oo and 
the population becomes extinct. 

Although the above model is a simplified 
one, it still reflects the actual situation in 
some cases. However, practically all mo- 
dels describing real phenomena and proces- 
ses are nonlinear, and instead of the differ- 
ential equation (32) we are forced to con- 
sider an equation of the type 


dx 
at | (x), 
where f (xz) is a nonlinear function, say the 


equation 


dz 
aol (x) = at — bxz?, 


where 

a>0,b> 0. 

Assuming that x = z, at ¢ = ¢, and solving 
the last equation, we get 


_ x,a/b 
*O= Zr @p—mewt—aeaan 8) 


We see that as t — oo the number of indiv- 
iduals in the population, z (t), tends to 
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Fig. 6 


a/b. Two cases are possible here: a/b > zy 
and a/b <1 zy. The difference between the 
two is clearly seen in Figure 6. Note that 
formula (33) describes, for one thing, the 
populations of fruit pests and some types 
of bacteria. 

If we consider several coexisting species, 
say, big and small fish, where the small fish 
serve as prey for the big, then by setting up 
differential equations for each species we 
arrive at a system of equations 


dz; 
Goa hi (yt), t= 1, 2% 


Let us study in greater detail the two-spe- 
cies predator vs. prey model, first intro- 
duced by the Italian mathematician Vito 
Volterra (1860-1940) to explain the oscilla- 
tions in catch volumes in the Adriatic Sea, 
which have the same period but differ in 
phase. 
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Let x be the number of big fish (predators) 
that feed on small fish (prey), whose num- 
ber we denote by y. Then the number of 
predator fish will grow as long as they have 
sufficient food, that is, prey fish. Finally, a 
situation will emerge in which there will 
not be enough food and the number of big 
fish will diminish. As a result, starting from 
a certain moment the number of small fish 
begins to increase. This, in turn, assists a 
new growth in the number of the big spe- 
cies, and the cycle is repeated. The model 
constructed by Volterra has the form 


rm —=cr— dy, (39) 


where a, b, c, and d are positive constants. 

In Eq. (34) for the big fish the term bzy 
reflects the dependence of the increase in 
the number of big fish on the number of 
small fish, while in Eq. (35) the term 
—dzxy expresses the decrease in the number 
of small fish as a function of the number of 
big fish. 

To make the study of these two equations 
more convenient we introduce dimension- 
less variables: 


u(t)=“2, v(t)=—y, t=ct, a=ale. 
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As a result the differential equations (34) 
and (35) assume the form 

u’ = au (v — 1), vb =vi(i — v), (36) 
with @ positive and the prime standing for 
differentiation with respect to T. 

Let us assume that at time t = ty, the 

number of individuals of both species is 
known, that is, 
ATG) Sy VG) SV (37) 
Note that we are interested only in positive 
solutions. Let us establish the relation- 
ship between u and v. To this end we divide 
the first equation in system (36) by the sec- 
ond and integrate the resulting differential 
equation. We get 


av + u — In vu = avg + Uy — In VU, 
af 


where H is a constant determined by the 
initial conditions (387) and parameter a. 

Figure 7 depicts the dependence of u 
on v for different values of H. We see that 
the (u, v)-plane contains only closed curves. 
Let us now assume that the initial val- 
ues U, and vy, are specified by point A 
on the trajectory that corresponds to the 
value H = H,. Since u, >1 and vy <1, 
the first equation in (36) shows that at first 
variable u decreases. The same is true of 
variable v. Then, as u approaches a value 
close to unity, v’ vanishes, which is fol- 
lowed by a prolonged period of time t in the 
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Fig. 8 


course of which variable v increases. When 
v becomes equal to unity, uw’ vanishes and 
variable uw begins to increase. Thus, both 
u and v traverse a closed trajectory. This 
means that the solutions are functions peri- 
odic in time. The variable u does not achieve 
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its maximum when v does, that is, oscil- 
lations in the populations occur in differ- 
ent phases. A typical graph of u and v as 
functions of time is shown in Figure 8 (for 
the case where v,>>1 and uy < 1). 

In conclusion we note that the study of 
communities that interact in a more com- 
plex manner provides results that are more 
interesting from the practical standpoint 
than those obtained above. For example, if 
two populations fight for the same source of 
food (the third population), it can be demon- 
strated that one of them will become 
extinct. Clearly, if this is the third popula- 
tion (the source of food), then the other two 
will also become extinct. 


1.9 A Problem from the 
Mathematical Theory of Epidemics 


Let us consider a differential model encoun- 
tered in the theory of epidemics. Suppose 
that a population consisting of N individ- 
uals is split into three groups. The first in- 
clude individuals that are susceptible to a 
certain disease but are healthy. We denote 
the number of such individuals at time 
t by S (t). The second group incorporates 
individuals that are infected, that is, they 
are ill and serve as a source of infection. 
We denote the number of such idividuals 
in the population at time ¢ by J (¢). Final- 
ly, the third group consists of individuals 
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who are healthy and immune to the di- 
sease. We denote the number of such indivi- 
duals at time ¢ by R (t). Thus, 


SQ +10 +RO)=N. (38) 


Let us also assume that when the number 
of infected individuals exceeds a certain 
fixed number /*, the rate of change in the 
number of individuals susceptible to the 
disease is proportional to the number of 
such individuals. As for the rate of change 
in the number of infected individuals that 
eventually recover, we assume that it is 
proportional to the number of infected in- 
dividuals. Clearly, these assumptions sim- 
plify matters and in a number of cases they 
reflect the real situation. Because of the 
first assumption, we suppose that when the 
number of infected individuals J (t) is 
greater than /*, they can infect the individ- 
uals susceptible to the disease. This means 
that we have allowed for the fact that the 
infected individuals have been isolated for 
a certain time interval (as a result of quar- 
antine or because they have been far from 
individuals susceptible to the disease). 
We, therefore, arrive at the following dif- 
ferential equation 


ds = —aS if I(t) > /*, 
dt. t Q, if, TAS i 4, 


Now, since each individual susceptible 
to the disease eventually falls ill and be- 


(39) 
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comes a source of infection, the rate of 
change in the number of infected individuals 
is the difference, per unit time, between the 
newly infected individuals and those that 
are getting well. Hence, 


dl aS—Bl if J (t) > I*, a 
dt — fl if I(t) < I*. (40) 
We will call the proportionality factors 
a and B the coefficients of illness and recovery. 
Finally, the rate of change of the number 
of recovering individuals is given by the 
equation dR/dt = Bl. 

For the solutions of the respective equa- 
tions to be unique, we must fix the initial 
conditions. For the sake of simplicity we 
assume that at time ¢t = 0 no individuals 
in the population are immune to the disease, 
that is, R (0) = 0, and that initially the 
number of infected individuals was J (0). 
Next we assume that the illness and recov- 
ery coefficients are equal, that is, a = 
(the reader is advised to study the case 
where a = f). Hence, we find ourselves in 
a situation in which we must consider two 
cases. 

Case 1. I (O)< J*. With the passage of 
time the individuals in the population will 
not become infected because in this case 
dS/dt = 0, and, hence, in accordance with 
Eq. (38) and the condition R (0) = 0, 
we have an equation valid for all values 


= 
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Fig. 9 


of t: 
S(t) = § (0) = N —TI (0). 


The case considered here corresponds to the 
situation when a fairly large number of in- 
fected individuals are placed in quaran- 
tine. Equation (40) then leads us to the fol- 
lowing differential equation 


di 

a7 —aal. 

This means that / (t) = J (0) e-* and, 
hence, 


R(t) = N — S (t) —T (@) 

= [ (0) [14 — e-], 

Figure 9 provides diagrams that illustrate 
the changes with the passage of time in the 
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number of individuals in each of the three 
groups. 

Case 2. I (0) > I*. In this case there must 
exist a time interval O< t< 7 in which 
I (t) > JI* for all values of ¢, since by its 
very meaning the function J = J (é) is 
continuous. This implies that for all ts 
belonging to the interval [0, 7] the disease 
spreads to the individuals susceptible to it. 
Equation (39), therefore, implies that 


S (t) = S (je 
for O< t< T. Substituting this into Eq. 
(40), we arrive at the differential equation 


Stal =a (0) e-. (1) 


If we now multiply both sides of Eq. (41) 
by e@f, we get the equation 


(let) =a8 (0). 


Hence, Je*# = aS (O)t +C and, there- 
fore, the set of all solutions to Eq. (41) is 
given by the formula 


I (t) = Ce-*t + aS (0) te-. (42) 


Assuming that ¢ = 0, we get C = IJ (QO) and, 
hence, Eq. (42) takes the form 
I (t) = (I (0) + aS (0) é] e-@, (43) 


with O< t< T. 
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We devote our further investigations to 
the problem of finding the specific value of 
T and the moment of time ¢,,,; at which the 
number of infected individuals proves 
maximal. 

The answer to the first question is impor- 
tant because starting from 7 the individuals 
susceptible to the disease cease to become 
infected. If we turn to Eq. (43), the afore- 
said implies that at t = T its right-hand 
side assumes the value J*, that is, 


I* = {I (0) + aS (0) Tle-@?, (44) 
But 
S (T) = lim S (t) = S (0) 

too 
is the number of individuals that are suscep- 
tible to the disease and yet avoid falling 


ill; for such individuals the following chain 
of equalities holds true: 


S (T) = S (co) = S (0) e-F, 
Hence, 


5) 
T= In F (cay * (45) 


Thus, if we can point to a definite value 
of S (co), we can use Eq. (45) to predict the 
time when the epidemic will stop. Substitut- 
ing T given by (45) into Eq. (44), we arrive 
at the equation 


S (0 ore) 
= [10)+5(Q mo |S , 


4—0770 
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or 
I* I (0) 


7 5 (0) 
Sto) SO) tS) 


which can be rewritten in the form 


i* __ £0) 


Since [* and all the terms on the right-hand 
side of (46) are known, we can use this equa- 
tion to determine S (oo). 

To answer the second question, we turn 
to Eq. (43). In accordance with the posed 
question, this equation yields 


=las (0) —al (0) — a2S (0) t] e-%# = 0. 


The moment in time at which / attains its 
maximal value is given by the following for- 
mula: 

4 I (0) 
tmax = =| ~ (0) ° 
If we now substitute this value into Eq. (43), 
we get 


Imax = S (0) exp{—[1—Z (0)/S (0)} 
=8 (tmax)- 


This relationship shows, for one, that at 
time ty,x the number of individuals sus- 
ceptible to the disease coincides with the 
number of infected individuals. 
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Fig. 10 


But if ¢ > 7, the individuals susceptible 
to the disease cease to be infected and 


L (t)= [¥#p-a(t- tT). 


Figure 10 gives a rough sketch of the dia- 
grams that reflect the changes with the pas- 
sage of time of the number of individuals 
in each of the three groups considered. 


1.10 The Pursuit Curve 


Let us examine an example in which differ- 
ential equations are used to choose a cor- 
rect strategy in pursuit problems. 


4&* 
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0 


Fig. 41 


We assume that a destroyer is in pursuit 
of a submarine in a dense fog. At a certain 
moment in time the fog lifts and the sub- 
marine is spotted floating on the surface 
three miles from the destroyer. The destroy- 
er’s speed is twice the submarine’s speed. 
We wish to find the trajectory (the pursuit 
curve) that the destroyer must follow to 
pass exactly over the submarine if the latter 
immediately dives after being detected 
and proceeds at top speed and in a straight 
line in an unknown direction. 

To solve the problem we first introduce 
polar coordinates, r and 9, in such a manner 
that the pole O is located at the point at 
which the submarine was located when dis- 
covered and that the point at which the 
destroyer was located when the submarine 
was discovered lies on the polar axis r 
(Figure 11). Further reasoning is based on 
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the following considerations. First, the de- 
stroyer must take up a position in which it 
will be at the same distance from pole O 
as the submarine. Then it must move about 
O in such a way that both moving objects 
remain at the same distance from point O all 
the time. Only in this case will the destroyer 
eventually pass over the submarine when 
circling pole O. The aforesaid implies that 
the destroyer must first go straight to point 
O until it finds itself at the same distance 
x from O as the submarine. 

Obviously, distance x can be found either 
from the equation 
x 3— 2 
‘vv 
or from the equation 
x 342 


‘vy ee 2 

where v is the submarine’s speed and 2v 
the destroyer’s speed. Solving these equa- 
tions, we find that this distance iseither one 
mile or three miles. 

Now, if the two still have not met, the 
destroyer must circle pole O (clockwise or 
counterclockwise) and head away from the 
pole at a speed equal to that of the subma- 
rine, v. Let us decompose the destroyer’s ve- 
locity vector (of length 2v) into two compo- 
nents, the radial component v, and the tan- 
gential component v, (Figure 11). 
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The radial component is the speed at 
which the destroyer moves away from pole 
O, that is 


dr 
Y, = ad’ 
while the tangential component is the li- 
near speed at which the destroyer circles 
the pole. The latter component, as is known, 
is equal to the product of the angular ve- 
locity d0@/dé and radius r, that is, 


v Baste AU 
t~" dt ° 


But since v, must be equal to v, we have 
v= V (vy? — v2 = V 3v. 


Thus, the solution of the initial problem 
is reduced to the solution of a system of two 
differential equations, 


dr dé ry 
am Ta HV 8p, 


which, in turn, can be reduced to a single 
equation, dr/r = d0/ 3, by excluding the 
variable ¢. 


Solving the last differential equation, we 
find that 


r=CelV 3, 


with C an arbitrary constant. 
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If we now allow for the fact that destroyer 
beginsits motion about pole O starting from 
the polar axis r at a distance of zx miles away 
from O, thatis,r = 1 at®0 = O andr = 3 at 
§ =—az, we conclude that in the first case 
C = 1 and in the second C = 3e*/V3, Thus, 
to fulfill its mission, the destroyer must 
move two or six miles along a straight line 
toward the point where the submarine was 
discovered and then move in the spiral 


r=e9V3 or the spiral r = 3e0+™/V3, 


1.141 Combat Models 


During the First World War the British 
engineer and mathematician Frederick 
W. Lanchester (1868-1946) constructed sev- 
eral mathematical models of air battles. 
Later these models were generalized so as to 
describe battles involving regular troops or 
guerilla forces or the two simultaneously. 
Below we consider these three models. 
Suppose that two opposing forces, xz 
and y, are in combat. We denote the num- 
ber of personnel at time ¢ measured in days 
starting from the first day of combat opera- 
tions by z (t) and y (t), respectively. It is 
the number of personnel that will play a 
decisive part in the construction of these 
models, since it is difficult, practically 
speaking, to specify the criteria that would 
allow taking into account, when comparing 
the opposing sides, not only the number of 
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personnel but also combat readiness, level 
of military equipment, preparedness and 
experience of the officers, morale, and a 
great many other factors. 

We also assume that x(t) and y(t) change 
continuously and, more than that, are 
differentiable as functions of time. Of course, 
these assumptions simplify the real sit- 
uation because z (t) and y (¢) are integers. 
But at the same time it is clear that if the 
number of personnel on each side is great, 
an increase by one or two persons will have 
an infinitesimal effect from the practical 
standpoint if compared to the effect pro- 
duced by the entire force. Therefore, we may 
assume that during smal] time intervals 
the change in the number of personnel is 
also small (and does not constitute an inte- 
gral number). The above reasoning is, of 
course, insufficient to specify concrete for- 
mulas for x (¢) and y (¢) as functions of t, but 
we can already point to a number of factors 
that enable describing the rate of change 
in the number of personnel on the two 
sides. Precisely, we denote by OLR the rate 
at which side z suffers losses from disease 
and other factors not related directly to 
combat operations, and by CLR the rate 
at which it suffers losses directly in combat 
operations involving side y. Finally by RR 
we denote the rate at which reinforcements 
are supplied to side z. It is then clear that 
the total rate of change of x (¢) is given by 
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the equation 


dz(t) 
dt 


A similar equation holds true for y (2). 
The problem now is to find the appropriate 
formulas for the quantities OLR, CLR, 
and RR and then examine the resulting 
differential equations. The results should 
indicate the probable winner. 

We use the following notation: a, b, 
c, d, g, and h are nonnegative constants 
characterizing the effect of various factors 
on personnel losses on both sides (x and y), 
P (t) and Q (t) are terms that allow for the 
possibility of reinforcements being supplied 
to x and y in the course of one day, and 
Xo and y, are the personnel in x and y prior 
to combat. We are now ready to set up the 
three combat models suggested by Lan- 
chester.* The first refers to combat opera- 
tions involving regular troops: 


te () — _ ax (t)—by (t) +P (2), 


su) —cex (f)—dy (1) + Q(t). 


—(OLR+CLR)+RR. (47) 


* Examples of combat operations are given in the 
article by C.S. Coleman, “Combat models” (see 
Differential Equation Models, eds. M. Braun, 
C.S. Coleman, and D.A. Drew, New York: Sprin- 
ger, 1983, pp. 109-131). The article shows how 
the examples agree with the models considered, 
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In what follows we call this system the A- 
type differential system (or simply the A- 
type system). 

The second model, specified by the equa- 
tions 


2) — — ax (t)—ga (t) y(t) +P (2), 
“HO — _ dy (t)—ha(t) y(t) + Q(t) 


describes combat operations involving only 
guerilla forces. We will call this system the 
B-type system. Finally, the third model, 
which we call the C-type system, has the 
form 


at) — — ax (t)— ga (t)y (t) + P (2), 
oo -- cx (t)— dy (t) +Q (t) 


and describes combat operations involving 
both regular troops and guerilla forces. 

Each of these differential equations 
describes the rate of change in the number 
of personnel on the opposing sides as a func- 
tion of various factors and has the form 
(47). Losses in personnel that are not direct- 
ly related to combat operations and are de- 
termined by the terms —az(t) and 
—dy (t) make it possible to describe the 
constant fractional loss rates (in the absence 
of combat operations and reinforcements) 
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via the equations 


1 dz 4 dy 


ga Oo pap 

If the Lanchester models contain only 
terms corresponding to reinforcements and 
losses not associated with combat opera- 
tions, this means that there are no com- 
bat operations, while the presence of the 
terms —by (t), —czx (t), —gzx (t) y (t), and 
—hz (t) y (t) means that combat operations 
take place. 

In considering the A-type system, we 
assume, first, that each side is within the 
range of the fire weapons of the other and, 
second, that only the personnel directly 
involved in combat come under fire. Under 
such assumptions Lanchester introduced the 
term —by (¢) for the regular troops of side 
x to reflect combat losses. The factor b 
characterizes the effectiveness of side y in 
combat. Thus, the equation 


1 dz 

7a 

shows that constant b measures the average 
effectiveness of each man on side y. The 
same interpretation can be given to the term 
—czx (t). Of course, there is no simple way 
in which we can calculate the effectiveness 
coefficients b and c. One way is to write 
these coefficients in the form 


b= PyPy: C = PlyPx, (48) 
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where r, andr, are the coefficients of fire- 
power of sides y and xz, respectively, and p, 
and p, are the probabilities that each shot 
fired by sides y and x, respectively, proves 
accurate. 

We note further that the terms correspond- 
ing to combat losses in the A-type system 
are linear, whereas in the B-type system 
they are nonlinear. The explanation lies in 
the following. Suppose that the guerilla 
forces amounting to x (t) personnel occupy a 
certain territory R and remain undetected 
by the opposing side. And although the lat- 
ter controls this territory by firepower, it 
cannot know the effectiveness of its ac- 
tions. It is also highly probable that the 
losses suffered by the guerilla forces x are, 
on the one hand, proportional to the num- 
ber of personnel z(t) on R and, on the 
other, to the number y (t) of personnel of 
the opposing side. Hence, the term corre- 
sponding to the losses suffered by the gueril- 
la forces x has the form —gz (t) y (¢), where 
the coefficient reflecting the effectiveness of 
combat operations of side y is, in general, 
more difficult to estimate than the coef- 
fiient b in the first relationship in (48). 
Nevertheless, to find g we can use the fire- 
power coefficient r, and also allow for the 
ideas expounded by Lanchester, according 
to which the probability of an accurate 
shot fired by side y is directly proportional 
to the so-called territorial effectiveness A,, 
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of a single shot fired by y and inversely 
proportional to the area A, of the territory 
R occupied by side xz Here by A,, we 
denote the area occupied by a single gueril- 
la fighter. Thus, with a high probability 
we can assume that the formulas for 
finding g and h are 


ems Ary a= Ap x 
§&=lry, Ax ’ h=r, Ay ° (49) 


Below we discuss in greater detail each of 
the three differential models. 

Case A (differential systems of the A-type 
and the quadratic law). Let us assume that 
the regular troops of the two opposing 
sides are in combat in the simple situation 
in which the losses not associated directly 
with combat operations are nil. If, in addi- 
tion, neither side receives any reinforce- 
ment, the mathematical model is reduced 
to the following differential system: 


Laan OW 
aay Gp ee. (50) 


Dividing the second equation by the first, 
we find that 


dy cz 
=F. (51) 


Integrating, we arrive at the following re- 
lationship: 


b {y? ¢) — y3] = e [z? (t) — 2). (52) 
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AG K?70, y wins 
K =0, even 


VX /b KO & wens 


V-“/4e T(t) 
Fig. 42 


This relationship explains why system (50) 
corresponds to a model with a quadratic 
law. If by K we denote the constant quanti- 
ty by} — cxj, the equation 


by? — cr? = K (53) 


obtained from Eq. (52) specifies a hyperbola 
(or a pair of straight lines if K = 0) and 
we can Classify system (50) with great pre- 
cision. We can call it a differential system 
with a hyperbolic law. 

Figure 12 depicts the hyperbolas for differ- 
ent values of K; for obvious reasons we 
consider only the first quadrant (« > 0, 
y > 0). The arrows on the curves point in 
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the direction in which the number of per- 
sonnel changes with passage of time. 

To answer the question of who wins in 
the constructed model (50), we agree to say 
that side y (or x) wins if it is the first to 
wipe out the other side z (or y). For exam- 
ple, in our case the winner is side y if 
K > 0 since in accordance with Eq. (53) 
the variable y can never vanish, while the 
variable x vanishes at y (t) = VY K/b. Thus, 
for y to win a victory, it must strive for a 
situation in which K is positive, that is, 


by? > cz. (54) 

Using (48), we can rewrite (54) in the form 
Yo \?, tx Px 

Cs ae Py ~ (59) 


The left-hand side of (55) demonstrates that 
changes in the personnel ratio yo/z» give 
an advantage to one side, in accordance 
with the quadratic law. For example, a 
change in the y)/z, from one to two gives y 
a four-fold advantage. Note also that 
Eq. (53) denotes the relation between the 
personnel numbers of the opposing sides but 
does not depend explicitly on time. To de- 
rive formulas that would provide us with a 
time dependence, we proceed as follows. 
We differentiate the first equation in (50) 
with respect to time and then use the second 
equation in this system. As a result we 
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Vi/o 


0 


Fig. 13 


arrive at the following differential equation: 


St —bex =0. (56) 


If for the initial conditions we take 


dx 

x (0) = 2p, leo — bY, 
then the solution to Eq. (56) can be ob- 
tained in the form 
x (t) = x) cos Bt — yy, sin Bi, (57) 
where 8 = Vbe and y = Vodice. 

In a similar manner we can show that 
y (t) = y, cos Bt — a sin pe. (98) 


Figure 13 depicts the graphs of the func- 
tions specified by Eqs. (57) and (58) in 
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the case where K >0O (i.e. by} > cx?, or 
Yo > Zo). 

We note in conclusion that for side y to 
win a victory it is not necessary that y, be 
greater than zx). The only requirement is 
that yy, be greater than 7. 

Case B (differential systems of the B-type 
and the linear law). The dynamical equa- 
tions that model combat operations between 
two opposing sides can easily be solved, 
as in the previous case, if we exclude the 
possibility of losses not associated with 
combat and if neither side receives rein- 
forcements. Under these limitations the 
B-type differential system assumes the form 


dz dy 
Gp Bt) Gp = hay. (99) 


Dividing the second equation in (59) by 
the first, we get the equation 


dy he 

dx g’ 

which after being integrated yields 

gly (t) — yo) = h Iz (t) — ol. (60) 


This linear relationship explains why the 
nonlinear system (59) corresponds to a mo- 
del with a linear law for conducting combat 
operations. Equation (60) can be rewritten 
in the form 


5—0770 
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yi) L>0, ywins 


L=0, ever 
LIg L<0, © Wins 
Q 7 
-L/h Z(t) 
Fig. 14 


with L = gy) — hx. This implies, for one, 
that if Z is positive, side y wins as a result 
of combat operations, while if Z is nega- 
tive, side x wins. 

Figure 14 provides a geometrical inter- 
pretation of the linear law (61) for different 
values of L. 

Let us now study in greater detail the 
situation where one of the sides wins. We 
assume that this is side y. Then, as we al- 
ready know, gy) — hz, must be positive, or 


Yo sh 

% 7 g° 

If we now turn to (49), we see that side y 
wins if 


Yo_ reApyAx 62 
Zo 7 ryApyAy ° we 
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Thus, the strategy of y is to make the ratio 
Yo/Xy as large as possible and the ratio 
A,/Ay as small as possible. From the practi- 
cal standpoint it iS more convenient to 
write inequality (62) in the form 


AyVo reAp x 
AxXo ryAry 


We see that in a certain sense the products 
A,yy, and A,xz, constitute critical values. 

We note, finally, that by combining 
(61) and (59) we can easily derive formulas 
that give the time dependence of changes 
in the personnel numbers of both sides. 

Case C (differential systems of the C-type 
and the parabolic law). In the C-model the 
guerilla forces face regular troops. We in- 
troduce the simplifying assumptions that 
neither is supplied with reinforcements and 
that the losses associated directly with com- 
bat operations are nil. In this case we have 
the following system of differential equa- 
tions: 


dz dy 
ar — §xry, ‘dt = — CL, (63) 


where x (t) is the number of personnel in 
the guerilla forces, and y (t) the number of 
personnel in the regular troops. Dividing 
the second equation in (63) by the first, we 
arrive at the differential equation 


dy € 
dx 
5* 
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YE)A 10, ¥ wins M=0 even 


M<O x2 wins 


Vay 


-M/2c Z(t) 


Fig. 15 


Integrating with the appropriate limits 
yields the following relationship: 


gy? (t) == 2er (t) + M, (64) 


where M = gy? — 2czy. Thus, the system 
of differential equations (63) corresponds to 
a model with a parabolic law for conducting 
combat operations. The guerilla forces win 
if M is negative; they are defeated if VW is 
positive. 

Figure 15 depicts the parabolas defined 
via Eq. (64) for different values of M. 
Experience shows that regular troops can 
defeat guerilla forces only if the ratio yo/zx, 
considerably exceeds unity. Basing our rea- 
soning on the parabolic law for conducting 
combat operations and assuming ™ posi- 
tive, we conclude that the victory of regular 
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troops is guaranteed if (y,/xz,)* is greater 
than (2c/g) 2;'. If we allow for (48) and 
(49), we can rewrite this condition in the 
form 


(2)? ote Aspe 4 
Ly Ty Ary te 


1.12 Why Are Pendulum Clocks 
Inaccurate? 


To answer this question, let us consder an 
idealized model of a pendulum clock con- 
sisting of a rod of length J and a load of mass 
m attached to the lower end of the rod 
(the rod’s mass is assumed so small that it 
can be ignored in comparison to m) (Fig- 
ure 16). Ifthe load is deflected by an angle a 
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and then released, in accordance with 
energy conservation we obtain 


mv2 
2 


= mg (lcos8—I cosa), (65) 


where v is the speed at which the load 

moves, and gis the acceleration of gravity. 
Since s=10, we have v= ds/dt = 

1 (d8/dt) and (65) leads us to the following 

differential equation: 

4 ( 2)’ =¢ (cos 0— cosa). (66) 

If we now allow for the fact that 6 decreases 


with the passage of time ¢ (for small t's), 
we can rewrite Eq. (66) in the form 


_ aes dé 
28 = cos0—cosa 


If we denote the period of the pendulum 
by 7, we have 


0 
i ee | VSS 
4 2e J Ycos0—cosa ” 
or 
Te do 
28 d Y cos 8 —cos a Sa 


The last formula shows that the period of 
the pendulum depends on a. This is the 
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main reason why pendulum clocks are inac- 
curate since, practically speaking, every 
time that the load is deflected to the ex- 
treme position the deflection angle differs 
from a. 

We note that formula (67) can be written 
in a simpler form. Indeed, since 


cos0=1—2Qsin?s, cosa=1—2 sin? > : 


we have 


—_— & 
= a ee 
' vs Verte 
=2f +t | panty 


A wate 2 k2?—sin? — 


with k = sin (a@/2). 

Now instead of variable 6 we introduce 
a new variable gq via the formula 
sin (0/2) = k sin g. This implies that when 
@ increases from 0 to a, the variable 
grows from 0 to 2/2, with 


(68) 


4 ) 
F COS d§ =k cos gdgq, 
or 


2k cos p dQ _ 2 2 sin ‘i —sint 
aaa, ada  Vt— sin? @ 
COS -5 ? 
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The last relationship enables us to rewrite 
formula (68) in the form 


1/2 


od dq 
ra4y/ 2 


V 1—? sin? 


“T 
= V =F, /2), 


where the function 


w 
_ dp 
F (, ») = V1—k? sin? — hk? sin? rs 


is known as the elliptic integral of the first 
kind, differing from the elliptic integral of 
the second kind 


p 
E(k, 1) = \ Vi-sin?@ dg. 


0 


Elliptic integrals cannot be expressed in 
terms of elementary functions. Therefore, 
all further discussion of the pendulum prob- 
lem will be related to an approach consid- 
ered when conservative systems are studied 
in mechanics. Here we only note that the 
starting point in our studies will be the 
differential equation 
Se tksind=0, k= Veil, 
which can be obtained from Eq. (60) by 
differentiating with respect to 7. 
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1.13 The Cycloidal Clock 


We have established that clocks with ordi- 
nary (circular) pendulums are inaccurate. 
Is there any pendulum whose period is in- 
dependert of the swing? This problem was 
first formulated and solved as early as the 
17th century. Below we give its solution, 
but first let us turn to the derivation of the 
equation of a remarkable curve, which 
Galileo Galilei called the cycloid (from the 
Greek word for “circular’). This is a plane 
curve generated by a point on the circum- 
ference of a circle (called the generating 
circle) as it rolls along a straight line with- 
out slippage. 

Suppose that the x axis is the straight 
line along which the generating circle rolls 
and that the radius of this circle is r (Fig- 
ure 17). Let us also assume that initially 
the point that traces the cycloid is at the 
origin and that after the circle turns through 
an angle 8, the point occupies position M. 
Then, on the basis of geometrical reasoning 
we obtain 


z= OS = OP — SP, 
y = MS = CP —CN. 
But 


OP = MP =r, SP = MN =rsin 89, 
CP =r, CN = rcos 0, 
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Fig. 17 


Hence, parametrically the cycloid is speci- 
fied by the following equations: 


x=r(§@—sin®8@), y=r(t — cos 86). (69) 


If we exclude parameter @ from these equa- 
tions, we arrive at the following equation 
of a cycloid in a rectangular Cartesian co- 
ordinate system Ozzy: 


x =r cos"! | —# )- V 2ry—y?. 


The very method of constructing a cy- 
cloid implies that the cycloid consists of con- 
gruent arcs each of which corresponds to a 
full revolution of the generating circle.* 


* The reader may find many interesting facts about 
the cycloid and related curves in G.N. Berman’s 
book The Cycloid (Moscow: Nauka, 1980) (in 
Russian). 
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The separate arcs are linked at points 
where they have a common vertical tan- 
gent. These points, known as cusps, corre- 
spond to the lowest possible positions occu- 
pied by the point on the generating circle 
that describes the cycloid. The highest 
possible positions occupied by the same 
point lie exactly in the middle of each arc 
and are known as the vertices. The distance 
along the straight line between successive 
cusps is 2nr, and the segment of that 
straight line between successive cusps is 
known as the base of an arc of the cycloid. 

The cycloid possesses the following prop- 
erties: 

(a) the area bounded by an arc of a cycloid 
and the respective base is thrice the area of 
the generating circle (Galileo's theorem); 

(b) the length of one cycloid arc is four 
times the length of the diameter of the gen- 
erating circle (Wren’s theorem). 

The last result is quite unexpected, since 
to calculate the length of such a simple 
curve as the circumference of a circle it is 
necessary to introduce the irrational number 
ma, whose calculation is not very simple, 
while the length of an arc of a cycloid is 
expressed as an integral multiple of the ra- 
dius (or diameter) of the generating circle. 
The cycloid has many other interesting 
properties, which have proved to be extreme- 
ly important for physics and engineering. 
For example, the profile of pinion teeth 
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Fig. 18 


and of many types of eccentrics, cams, and 
other mechanical parts of devices have the 
shape of a cycloid. 

Let us now turn to the problem that 
enabled the Dutch physical scientist, astro- 
nomer, and mathematician Christian Huy- 
gens (1629-1695) to build an accurate clock 
in 1673. The problem consists of building 
in the vertical plane a curve for which the 
time of descent of a heavy particle sliding 
without friction to a fixed horizontal line 
is the same wherever the particle starts on 
the curve (it is assumed that initially, at 
time ¢ = f,, the particle is at rest). Huy- 
gens found that the cycloid possesses such 
an isochronous (from the Greek words for 
“equal” and “time”) or tautochronous 
property (from tautos for “identical”). 

From the practical standpoint the prob- 
lem can be solved in the following manner. 
Let us assume that a trough in the form of a 
cycloid is cut out of a piece of wood as 
shown in Figure 18. We take a small metal 
ball and send it rolling down the slope. 
Ignoring friction, let us try to determine 
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the time it takes the ball to reach the low- 
est possible point, K, starting from point 
M, say. 

Let x, and y, be the coordinates of the 
initial position of the ball, i.e. point M, 
and @, the corresponding value of para- 
meter @. When the ball reaches a point N (6), 
the distance of descent along the vertical 
will be hk, which in view of Eq. (69) can be 
found in the following manner: 


h=y—y,=r (i — cos 8) 
— r(1— cos 6,) = r (cos 0, — cos 8). 


We know that the speed of a falling object 
is given by the formula 


v= V Qgh, 


where g is the acceleration of gravity. In 
our case the last formula assumes the form 


v= V 2gr (cos @, — cos 8) 


On the other hand, since speed is the deriv- 
ative of distance s with respect to time 7, 
we arrive at the formula 


— V 2gr (cos 0, —cos 6) 

Since for a cycloid ds = 2r sin (8/2) d@, we 
can rewrite this formula (which is actually 
a differential equation) in the form 


2r sin (8/2) dO 


ie a a 
V 2er (cos 0) — cos 8) 
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Integrating this equation within approp- 
riate limits, we obtain 


m 
= ( 2r sin (6/2) dé 
V 2egr (cos 05 —cos 9) 


_ wa 
=—2y 2{ d cos (0/2) 
_ g or: as * i 
0 \/ cos? 3-08 > 


=m VY r/g. 


Thus, the time interval 7 in the course 
of which the ball rolls down from point M 
to point K is given by the formula 


T=nVrig, 


which shows that period 7 is independent 
of 8,, that is, the period does not depend 
on the initial position of the ball, MM. 
Obviously, two balls that begin their mo- 
tion simultaneously from points M and N 
will roll down and find themselves at point 
K at the same moment of time. 

Since we agreed to ignore friction, in its 
motion down the slope the ball passes point 
K and, by inertia, continues its motion 
up the slope to point M, lying at the same 
level as point M. Proceeding then in the 
opposite direction and traversing its path, 
the ball completes a full cycle. This con- 
stitutes the motion of a cycloidal pendulum 
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Fig. 49 


with a period of oscillations 
T, = 4nV rig. (70) 


A peculiar feature of the cycloidal pen- 
dulum, which sets it apart from the simple 
(circular) pendulum, is that its period does 
not depend on the amplitude (or swing). 

Let us now show how an ordinary circu- 
lar pendulum can be made to move in a 
tautochronous manner without resorting to 
trough and similar devices with consider- 
able friction. To this end it is sufficient 
to make a template (say, out of wood) con- 
sisting of two semiarcs of acycloid witha 
common cusp (Figure 19). The template 
is fixed in the vertical plane and a string 
with a ball is suspended from the cusp O. 
The length of the string must be twice the 
diameter of the generating circle of the 
cycloid. 

If the ball on the string is deflected to an 
arbitrary point M, it begins to swing back 
and forth with a period that is independ- 
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y ent of the position 
of point M. Even 
if due to friction 
and air drag the 
amplitude of the 
oscillations dimin- 
ishes, the period 
remains unchan- 
A ged. For a circu- 
8 lar pendulum, 
C £ which moves along 
an arc of a cir- 
Fig. 20 cumference, the iso- 
chronous _— proper- 
ty is satisfied approximately only for 
small amplitudes, when the arc differs 
little from the arc of a cycloid. 
By way of an example let us study the 


small oscillations that a pendulum executes 


along the arc AB of a cycloid (Figure 20). 
If these oscillations are very small, the 
effect of the guiding template is practically 
nil and the pendulum oscillates almost like 
an ordinary circular pendulum with a string 
(the “rod”) whose length is 4r. The path AB 
of a cycloidal pendulum will differ little 
from the path CE of a circular pendulum 
with a string length equal to 4r. Hence, the 
period of small oscillations of a circular 
pendulum with a string length 1 = 4r is 
practically the same as that of a cycloidal 
pendulum of the same length. 


L=4r 
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Now, if in (70) we put r = 1/4, the period 
of small oscillations of a circular pendu- 
lum can be expressed in terms of the length 
of the string: 

T = 2nV lig. 
This formula is derived (in a different way) 
in high school physics. 

In conclusion we note that the cycloid 
is the only curve for which a particle mov- 
ing along it performs isochronous oscilla- 
tions. 


1.44 The Brachistochrone Problem 


The problem concerning the brachistochrone 
(from the Greek words for “shortest” and 
“time”), a curve of fastest descent, was pro- 
posed by the Swiss mathematician John 
Bernoulli (1667-1748) in 1696 as a challenge 
to mathematicians and consists of the fol- 
lowing. 

Take two points A and B lying in a ver- 
tical plane but not on a single vertical line 
(Figure 21). Among the various curves pass- 
ing through these two points we must find 


A 


Fig. 24 
6—0770 


82 Differential Equations in Applications 


Fig. 22 


the one for which the time required for a 
particle to fall from point A to point B 
along the curve under the force of gravity is 
minimal. 

The problem was tackled by the best 
mathematicians. It was solved by John 
Bernoulli himself and also by Gottfried W. 
von Leibniz (1646-1716), Sir Isaac Newton 
(1642-1727), Guillaume L’Hospital (1661- 
1704), and Jakob Bernoulli (1654-1705). 
The problem is famous not only from the 
general scientific viewpoint but also for 
being the source of ideas in a completely 
new field of mathematics, the calculus of 
variations. 

Solution of the brachistochrone problem 
can be linked with that of another problem 
originating in optics. Let us turn to Fig- 
ure 22, in which aray of light is depicted as 
propagating from point A to point P with 


Ch. 4. Construction of Differential Models 85 


a velocity v, and then from point P to 
point B in a denser medium with a lower ve- 
locity v,. The total time 7 required for the 
ray to propagate from point A to point B 
can, obviously, be found from the following 
formula: 


is aaa Vere Se. 


If we assume that ce ray of light propa- 
gates from point A to point B along this path 
in the shortest possible time 7, the deriva- 
tive d7/dz must vanish. But then 


x _— c—2x 
vy Vara vg YORE (ea) 
or 
sing, sin @, 
v1 vg 


The last formula expresses Snell’s well- 
known law of refraction, which initially 
was discovered in experiments in the form 
sin @,/sin a = a, with a constant. 

The above assumption that light chooses 
a path from A to B that would take the 
shortest possible time to travel is known as 
Fermat's principle, or the principle of least 
time. The importance of this principle lies 
not only in the fact that it can be taken as 
a rational basis for deriving Snell's law but 
also, for one, in that it can be applied to 
finding the path of a ray of light in a me- 
dium of variable density, where generally 


6e@ 
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Fig. 23 


the light travels not along straight-line 
segments. 

For the sake of an example let us con- 
sider Figure 23, which depicts a ray of light 
propagating through a layered medium. In 
each layer the speed of light is constant, 
but it decreases when we pass from an up- 
per layer to a lower layer. The incident ray 
is refracted more and more strongly as it 
passes from layer to layer and moves closer 
and closer to the vertical line. Applying 
Snell’s law to the interfaces between the 
layers, we get 


SINn@, _SIN@, _—s- SING, ~—s SIN 
M% «=©6hmyts—<i<i<i‘<‘<“<‘<z ™St”*”*i‘(‘ yS*C 


Now let us assume that the layer thick- 
nesses decrease without limit while the num- 
ber of layers increases without limit. Then, 
in the limit, the velocity of light changes 
(decreases) continuously and we conclude 


= 
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Fig. 24 
(see Figure 24) that 
a, (71) 


Vv 
with a = const. A similar situation is ob- 
served (with certain reservations) when a 
ray of Sun light falls on Earth. As the ray 
travels through Earth’s atmosphere of in- 
creasing density, its velocity decreases and 
the ray bends. 

Let us go back to the brachistochrone 
problem. We introduce a system of coor- 
dinates in the vertical plane in the way 
shown in Figure 25. We imagine that a 
ball (like a ray of light propagating in me- 
dia) is capable of choosing the path of de- 
scent from point A to point B with the short- 
est possible time of descent. Then, as fol- 
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Fig. 25 


lows from the above reasoning, formula (71) 
is valid. 

The energy conservation principle implies 
that the speed gained by the ball at a 
given level depends only on the loss of po- 
tential energy as the ball reaches the level 
and not on the shape of the trajectory fol- 
lowed. This means that v= V 2gy. 

On the other hand, geometric construc- 
tion enables us to show that 


sina=cosf = ee eae eee 
secB = 11+ tan?B 
ees eee 
ViFy? | 
Combining the last two relationships with 
(71) yields 


y [t+ (y')"J=¢, (72) 


= 
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This constitutes the brachistochrone equa- 
tion. We wish to demonstrate that only a 
cycloid can represent a brachistochrone. 
Indeed, since y’ = dy/dz, by dividing the 
variables in Eq. (72) we arrive at the equa- 
tion 

- y 1/2 
dx = ( c= dy. 
We introduce a new variable @ via the 
following relationship: 


( ear == tan Q. 


Then 
y=Csin?, dy=2C sin gcos 9g dq, 


dz = tan m dy = 2C sin? p dg 
=C (1— cos 29) dg. 


Integration of the last equation yields 
x = (29 —sin 29) + Cy, 


where, in view of the initial conditions, 
x=y=0 at g=0 and C, = 0. Thus, 


rs s (2m — sin 29), 
y=C sin? p= 5 (1 — cos 29). 


Assuming that C/2 =r and 2p = 0, we 
arrive at the standard parametric equations 
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of the cycloid (69). The cycloid is, indeed, a 
remarkable curve: it is not only isochro- 
nous—it is brachistochronous. 


1.15 The Arithmetic Mean, the 
Geometric Mean, and the 
Associated Differential Equation 


Let us consider the following curious prob- 
lem first suggested by the famous German 
mathematician Carl F. Gauss (4777-1855). 

Let m, and ny be two arbitrary positive 
numbers (m, >> 7,). Out of m, and ny, we 
construct two new numbers m, and n, that 
are, respectively, the arithmetic mean and 
the geometric mean of m, and ny. In other 
words, we put 

My + No 


Ne eg 2 ed V mana; re 


Treating m, and n, in the same manner as 


m, and nyo, we put 


my +My 
2 


m, = , Ny V myn. 


Continuing this process indefinitely, we ar- 
rive at two sequences of real numbers, {m;} 
and {n,} (k = 0, 4, 2, .), which, as can 
easily be demonstrated, are convergent. We 
wish to know the difference of the limits of 
these two sequences. 

An elegant solution of this problem 
amounting to setting up asecond-order linear 
differential equation is given below. It be- 
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longs to the German mathematician Carl 
W. Borchardt (4817-1880). Let a be the 
sought difference. [t obviously depends on 
m, and ny, a fact expressed analytically as 
follows: a =f (my, mo), Where f is a func- 
tion. The definition of a also implies that 
a=f(m,, n,). Now, if we multiply m, 
and n, by the same number k, each of the 
numbers m,, 1, Ms, No, introduced 
above, including a, will be multiplied by 
k. This means that a is a first-order homo- 
geneous function in m, and ny, and, hence, 


a = Mof (A, ro/myo) = mf (A, Ny/m). 


Introducing the notation x = no/my, 2 


n/m, *9 Yy — 1/f (4, No/M,), Wy —= 
4/f (4, n/m), we find that 
2 
Y= Ui Tae (73) 
Since z, is related to x via the equation 
2V a 
a {+2 ° 
we find that 
dz, 1—<2£ __ (ty — 23) (142)? 
dz (+422 Vx 2 (#2) 


On the other hand, Eq. (73) leads to the 
following relationship: 


dy 2 2 dy, dz, 
dz~ = (14-2)? Ww TLe dz, de ' 


= 
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Replacing dz,/dxz with its value given by 
the previous relationship and factoring out 
(x — x*), we find that 


d 2x (x—1 
(0) Fie 

d 
+(4+ 2) (a,— 24) Ge. 


Differentiating both sides of this equation 
with respect to x, we find that 


d | ( — x3) 54] 
dz 


x(x—1) 
al 44-2 J yes dy, dz, 
ae Essiad 


=2y1 1tz dz, dz 


3 dy 
+ (x, ~ 23) Es 


d 5 dz 
+(1-+2) 3| (a, —23) G2]. 


dz, 


By elementary transformation this equa- 
tion can be reduced to 


ae (@—2) #]—2y 


1—<z 


ge FON Ufo, ays 
~ (A+2) Va ass | ey 


dyy 
dz, 


]-eu.}. 


If we now replace z with z,, then z, will 
become z,. If we then replace x, with 2z,, 
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we find that z, is transformed into zz, etc. 
Hence, assuming that 


se | 2) 54] —y=a* (y), 


we arrive at the following formula: 


FM (eae Rc 
(A--2) Vx (4-+2y) Vay (144+22) Vtq 
m ae? (Ha). 


x ———_—_— 
(1-+ tp) V tp 


If we now send n to infinity, we find that 
1 — 2x, tends to zero and, hence, 


a* (y) = 0. 


This means that y satisfies the differential 
equation 

d2y 
dx? 
If we note that 


(x — x) +-(4 — 3x?) Sty =Q. (74) 


f? (Mo, No) 


a= f (IM, Mo) =Y Mg 


we can easily find the value of this num- 
ber. Indeed, since y must be the constant 
solution to Eq. (74), we find that only 
y =Oissuch a solution. Thus, the differ- 
ence of the limits of the sequences {m;,} and 
{n,} is zero. 

The differential equation (74) is remark- 
able not only because it enabled reducing 
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the initial problem to an obvious one but 
also because it is linked directly to the so- 
lution of the problem of calculating the pe- 
riod of small oscillations of a circular pen- 
dulum. 

As demonstrated earlier, the period of 
small oscillations of a circular pendulum 
can be found from the formula 


T=4VUgF (k, x/2), 


with 
/2 ig 
Ene) = Vishay 


It has been found that if O< k <1, then 


eee. eee 
V 1—k3 sin? p 


L ws) 


mfg) St 12XBEX ERX LX (221)? ) 
=F lt+2. 22 x 42 x 62X ... X (2n)? Ke) 
r= 


where 


— 42 B2X 52K... X (2N—1)2 on 
y=i+ > 22s 42 x 62 x X (2n)2 x? 
n=1 


is a solution to the differential equation (74). 
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1.16 On the Flight of an Object 
Thrown at an Angle to the Horizon 


Suppose that an object is thrown at an an- 
gle a to the horizon with an initial velocity 
Vy. We wish to derive the equation of the 
object’s motion that ignores forces of 
friction (air drag). 

Let us select the coordinate axes as shown 
in Figure 26. At an arbitrary point of the 
trajectory only the force of gravity P 
equal to mg, with m the mass of the object 
and g the acceleration of gravity, acts on 
the object. Hence, in accordance with 
Newton’s second law, we can write the 
differential equations of the motion of the 
object as projected on the x and y axes as 
follows: 


dr d@y 
mae aad a ee 
Factoring out m yields 
d?z d@y 
ao aa =F (79) 
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The initial conditions imposed on the ob- 
ject’s motion are 


x=0, y=0, < —y, cosa, 

t 
5 (76) 
Fp =o sina at t=0. 
Integrating Eqs. (75) and allowing for the 
initial conditions (76), we find that the 
equations of the object's motion are 


x = (vy cosa) t, y = (Uy sin a) t — gt?/2. 
(77) 


A number of conclusions concerning the 
character of the object’s motion can be 
drawn from Eqs. (77). For example, we can 
find the time of the object’s flight up to the 
moment when the object hits Earth, the range 
of the flight, the maximum height that the 
object reaches in its flight, and the shape 
of the trajectory. 

The first problem can be solved by finding 
the value of time ¢ at which y = 0. The 
second equation in (77) implies that this 
happens when 


, gt 
t | % sina —£-|=0, 


that is, either ¢ = 0 or t¢t = (2v, sin a)/g. 
The second value provides the solution. 

The second problem can be solved by 
calculating the value of x at a value of ¢ 
equal to the time of flight. The first equa- 
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tion in (77) implies that the range of the 
flight is given by the formula 


(v9 COS @) (2v9 3iINn a) _—svg sin 2a 
g gk 


This implies, for one thing, that the range 
is greatest when 2a = 90°, or @ = 45° In 
this case the range is v2/g. 

The solution to the third problem can be 
obtained immediately by formulating the 
maximum condition for y. But this means 
that at the point where y is maximal the 
derivative dy/dt vanishes. Noting that 


d ; 
<= —gt+ysina, 
we arrive at the equation —gt -++ v, sina =0, 
which yields 


t = (Vv, sin a)/y. 


Substituting this value of ¢ into the second 
equation in (77), we find that the maximum 
height reached by the object is v? sin? a/2g. 

The solution to the fourth problem has 
already been found. Namely, the trajectory 
is represented by a parabola since Eqs. (77) 
represent parametrically a parabola, which 
in rectangular Cartesian coordinates can 
be written as follows: 


gz” 
2v2 


y= zx tana— sec? a. 
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1.17 Weightlessness 


The state of weightlessness (zero g) can be 
achieved in various ways, although it is 
associated (consciously or subconsciously) 
with the “floating” of astronauts in the 
cabin of a spacecraft. 

Let us consider the following problem. 
Suppose that a person of weight P is stand- 
ing in an elevator that is moving down- 
ward with an acceleration wo = ag, with 
O<a<c1 and g the acceleration of grav- 
ity. Let us determine the pressure that 
the person exerts on the cabin’s floor and 
the acceleration that the elevator must 
undergo so that this pressure will vanish. 

Two forces act on the person in the eleva- 
tor (Figure 27): the force of gravity P and 
the force Q that the floor exerts on the per- 


Q 


Fig. 27 
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son (equal numerically to the force of pres- 
sure of the person on the floor). The differ- 
ential equation of the person’s motion can 
be written in the form 


d2 


Since d?z/dt? = wo = ag and m = P/g, we 
can rewrite Eq. (78) as follows: 


O=P— m Se =P (1—a). (79) 


Since 0 << a <1, we conclude that 0 < P. 
Thus, the pressure that the person exerts on 
the floor of the cabin of an elevator mov- 
ing downward is determined by the force 


Q=P(1— a). 


On the other hand, when the elevator is mov- 
ing upward with an acceleration o = ag, 
O<a<i, the pressure that the person 
exerts on the floor of the cabin is deter- 
mined by the force Q = P (1 + a). Let us 
now establish at what acceleration the pres- 
sure vanishes. For this it is sufficient to put 
Q = 0 in (79). We conclude that in this 
case a = 1, that is, for Q to vanish the 
acceleration of the elevator must be equal 
to the acceleration of gravity. 

Thus, when the cabin is falling freely 
with an acceleration equal to g, the pres- 
sure that the person exerts on the floor is 
nil. It is this state that is called weightless- 
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ness. In it the various parts of a person's 
body exert no pressure on each other, so 
that the person experiences extraordinary 
sensations. In the state of weightlessness 
all points of an object experience the same 
acceleration. 

Of course, weightlessness is experienced 
not only during a free fall in an elevator. 
For illustration let us consider the following 
problem. 

What must be the speed of a spacecraft 
moving around the Earth as an artificial 
satellite for a person inside it to be in the 
state of weightlessness? 

One assumption in this problem is that 
the spacecraft follows a circular orbit of 
radius r -+- h, where r is Earth’s radius, and 
h is the altitude at which the spacecraft 
travels (reckoned from Earth’s surface). The 
previous problem implies that in the state 
of weightlessness the pressure on the walls 
of the spacecraft is zero, so that the force Q 
acting on an object inside the spacecraft 
is zero, too. Hence, Q = 0. Let us now 
turn our attention to Figure 28. The x axis 
is directed along the principal normal n 
to the circular trajectory of the spacecraft. 
We use the differential equation of the mo- 
tion of a particle as projected on the prin- 
cipal normal: 


nr 
mov 
0 = >, Fras 


k=14 
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Fig. 28 


where 9 = r-+A, By Fin =F and F is 


directed along the rine pal normal to the 
trajectory. We get 


mv2 _P d2zx 


rh — Ne 
or the equation 

d2zx si 

dt? = rth 


Substituting this value of d?*z/di? into 
Eq. (78), we find that 


mv2 
=P 9. (80) 
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Here force P is equal to the force F of 
attraction to Earth, which, in accordance 
with Newton's law of gravitation, is inverse- 
ly proportional to the square of the dis- 
tance r-+h from the center of Earth, 
that is, 


where m is the mass of the spacecraft, and 
constant k can be determined from the 
following considerations. At Earth’s sur- 
face, where h = 0, the force of gravity F 
is equal to mg. The above formula then 
yields k = gr*. Hence, 


___ mgr? 
P= P= GLa? 
where g is the acceleration of gravity at 
Earth’s surface. 

If we now substitute the obtained value 
of P into (80) and note that Q = 0, we find 
that the required speed is given by the 
formula 


ae g 
very 


1.148 Kepler’s Laws of Planetary 
Motion 


In accordance with Newton’s law of gravi- 
tation, any two objects separated by a dis- 
tance r and having masses m and M are at- 


= 
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Fig. 29 


tracted with a force 
F= GmM (81) 


r2 9 


where G is the gravitational constant. 

Basing ourselves on this law, let us de- 
scribe the motion of the planets in the solar 
system assuming that m is the mass of a 
planet orbiting the Sun and M is the 
Sun’s mass. The effect of other planets on 
this motion will be ignored. 

Let the Sun be at the origin of the coor- 
dinate system depicted in Figure 29 and 
the planet be at time ¢ at a point with 
running coordinates zx and y. The attractive 
force F acting on the planet can be decom- 
posed into two components: one parallel! to 
the x axis and equal to # cosq and the 
other parallel to the y axis and equal to 
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F sin g. Using formula (81) and Newton’s 
second law, we arrive at the following 
equations: 


mz = — F cos Q=— one COS @, (82) 
my = —Fsing=— ome sin q. (83) 


Bearing in mind that sing = y/r and 
cos @ = az/r, we can rewrite Eqs. (82) and 
(83) as follows: 


kz si ky 

a Ths a acne a 

where constant k is equal to GM. 
Finally, allowing for the fact that r = 

VY x2 + y®, we arrive at the differential 

equations 


= kz ae ky 
OS Gap ye 9 4 aya 
(34) 


Without loss of generality we can assume 
that 


t=a,y=0,2=0,y=vn, att =0. (85) 


We see that the problem has been re- 
duced to the study of Eq. (84) under the 
initial conditions (85). The special fea- 
tures of Eqs. (84) suggest that the most con- 
venient coordinate system here is the one 
using polar coordinates: x =rcos@ and 
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y =rsing. Then 

x =Trcos 9 — (rsin 9) 0, 

y = r sin @ + (r cos @) 0, 

2 =rcos p — 2 (r sin @) © (86) 
— (rsin @) 9 — (r cos @) o?, 

y = rsin m+ 2 (r COS @) © 

+ (r cos q) © — (r sin g) @?. 

Hence, 

t= (r — rep?) cos mY — (2r @ 4. r@) sin q, 

y = (r — rg’) sin g + (279 + 1r@) cos 9g. 
Using the last two relationships, we can 


rewrite the differential equations (84) in 
the form 


(r~ ro?) cos ~ — (2re@ _ r@) sin @ 


k 
(r —rg’) sin @ -+ (2rp-++ r 0) cos © 
eiane ac Lea (88) 


Multiplying both sides of Eq. (87) by 
cos @, both sides of Eq. (88) by sin g, and 
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adding the products, we find that 


r— rg? = —kir*. (89) 
Multiplying both sides of Eq. (87) by 
sin g, both sides of Eq. (88) by cos g, and 
subtracting the second product from the 
first, we arrive at the equation 


2rp + rp = 0. (90) 
As for the initial conditions (85), in polar 
coordinates they assume the form 


r=a,9=0,r=0, o=v,/a at t=0. 
(91) 


Thus, we have reduced the problem of 
studying Eqs. (84) under the initial condi- 
tions (85) to that of studying Eqs. (89) and 
(90) under the initial conditions (91). We 
also note that Eq. (90) can be rewritten in 
the form 


d e 

—S (729) =0. (92) 
But Eq. (92) yields 

rp = Cy, (93) 


where C, is a constant possessing an in- 
teresting geometric interpretation. Precise- 
ly, suppose that an object moves from 


point P to point Q along the arc PO (Fig- 
ure 30). Let S be the area of the sector limit- 
ed by the segments OP and OQ and the 


= 
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Y 


0 


Fig. 30 


arc PQ. From the calculus course we know 
that 


@ 
4 
S => r2dq, 
0 
or dS = (1/2) r? dg. Hence, 
dS f dp  1.,° 
or ee a mee Pl (94) 


The derivative dS/dé constitutes what is 
known as the areal velocity, and since in 


view of (93) r’p is a constant, we conclude 
that the areal velocity is a constant, too. 
But this, in turn, means that the object 
moves in such a manner that the radius 
vector describes equal areas in equal time 
intervals. This law of areas constitutes one 
of the three Kepler laws. In full it can be 
formulated as follows: each planet moves 
along a plane curve around the Sun in such 
a manner that the radius vector connecting 
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the Sun with the planet describes equal areas 
in equal time intervals. 

To derive the next Kepler law, which 
deals with the shape of the planetary tra- 
jectories, we return to Eqs. (89) and (90) 
with the initial conditions (91) imposed on 
them. The initial conditions imply, for 


one, that r =a and O = voia at t = 0. 
But then condition (93) implies that C, = 
av,. Hence, 


rep = AVp, OF © =A ir". (95) 
This transforms Eq. (89) into 

"satu k 

BE age eee 


Assuming that r = p, we can rewrite this 
equation in the form 


dp dp drs dp ark 


‘dt dr dt = dye r3 “pe? 
or 


dp __ a*v§ k 
Par 


r3 r2 

Separating the variables in the last diffe- 
rential equation and integrating, we arrive 
at the oe relationship: 


p? a2y2 
oo r ae a3 C2. 


Since p = r=0 at r=a, we find that 


pe k 
C2= 3 a 
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Thus, we arrive at the equation 


r? k a2v2 i vz k 


—_—, 
— 


2 Tr 2 


or, if we consider only the positive value 
of the square root, 


dr 2k 2k a2y2 
a=V (%—--2)4+ 5 -S (96) 


Dividing Eq. (96) by Eq. (95), we find that 


dr |. EG ae 
de” V ar -+- 2Br—1, 
where 

1 2k __k 
O=  ady2? a a2v3 


The last equation can be integrated by sub- 
stituting 1/u for r. The result is 


a2v2/ke 


{+-ecos(p+Cs) ’ 


where e = Va + B?/B = av¥/k —1. The 
constant C, can be determined from the 
condition that r=a at » = 0. It is easy 
to verify that C, = 0. Thus, we finally have 

— atug/ie 

~ 4+.ecos @ (97) 
From analytic geometry we know that this 
is the equation of a conic in terms of polar 
coordinates, with e the eccentricity of the 
conic. The following cases are possible 


A 
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here: 
(1) an ellipse if e<4, or Ww < 2k/a, 
(2) a hyperbola if e>1, or v? > 2k/a, 
(3) a parabola if e = 1, or ve = 2k/a, 
(4) a circle ife = 0, or 2 = = hla. 


Astronomical observations have shown 
that for all the planets belonging to the 
solar system the value of v? is smaller 
han 2k/a. We, therefore, arrive at another 
of Kepler’s laws: the planets describe el- 
lipses with the Sun at one focus. 

Note that the orbits of the Moon and 
the artificial satellites of Earth are also 
ellipses, but in the majority of cases these 
ellipses are close to circles, that is, e differs 
little from zero. 

As for recurrent comets, like, say, Hal- 
ley’s comet, their orbits resemble “prolate” 
ellipses whose eccentricity is smaller than 
unity but very close to it. Say, Halley’s 
comet appears in Earth’s neighborhood 
approximately every 76 years. Its latest 
apparition was in the period between the 
end of 1985 and the beginning of 1986. 

Celestial bodies that move along para- 
bolic and hyperbolic orbits may be observed 
only once, since they never return to 
the same place. 

Let us now establish the physical mean- 
ing of eccentricity e. First we note that 


the components x and y of a planet’s vel- 
ocity vector V along the z and y axes, respec- 
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tively, and the size v of vector V satisfy 
the relationship 


v= x? Se y?, 
which when combined with (86) yields 
ites rp? che 7 


From this it follows that the kinetic energy 
of a planet of mass m is given by the for- 
mula 


my? => m (rg2-+r2), (98) 
The potential energy of a system is minus 
one times the amount of work needed to 
move the planet to infinity (where the po- 
tential energy is zero by convention). 
Hence, 


km km | km 
—-\ Fara =— (99) 


r r 


r 


If by E we denote the total energy of the 
system, which in view of energy conserva- 
tion must be constant, then formulas (98) 
and (99) yield 


5m [reg? tr] — A" =. (100) 


Assuming that » = 0 and combining (97) 
with (100), we get 


ue a2y2/k mrta®vg km ey, 
1te? 2r4 r : 
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Excluding r from the last two relationships, 
we find that 


> 2a%v2 
o- ices +E mk 


Equation (97) for the shape of the orbit 
finally assumes the form 


a2v2/k 
14- V1+E (2a2v2/mk?) cos 


This formula implies that the orbit is an 
ellipse, hyperboba, parabola, or circle if, 
respectively, E<0, E>0O, E =O, or 
E = —mk?/2a*v?. Thus, the shape of a 
planet’s orbit is completely determined by 
the value of E. Say, if we could impart 
such a “blow” to Earth that it would increase 
Earth’s total energy to a positive value, 
Earth would go over to a hyperbolic orbit 
and leave the solar system forever. 

Let us now turn to Kepler’s third law. 
This law deals with the period of revolu- 
tion of planets around the Sun. Taking 
into account the results obtained in deriv- 
ing the previous Kepler law, we naturally 
restrict our discussion to the case of ellip- 
tic orbits, whose equation in terms of Carte- 
sian coordinates, as is known, is 


x2 y? 
eta sh 


= 
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Fig. 34 
where (Figure 31) the eccentricity e = C/E 
with C? = &? — ?, so that 

e? == (6? — 1”)/6, 

or 

n’ = & (ft — e?). (101) 


Combining this with (97) and allowing for 
the properties of an ellipse, we conclude 


that 

4 f atugk avi/k \ atu 
Sy ( ite + 1—e )= ke (4 —e?) 
__ arugE? 
> ky? ? 
or 

2 

qe ie (102) 


k 
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Let us denote the period of revolution of a 
planet by 7. By definition, 7 is the time 
it takes the planet to complete one full or- 
bit about the Sun. Then, since the area lim- 
ited by an ellipse is m&n, we conclude, on 
the basis of (94) and (95), that mEn = 
av,T/2. Finally, taking into account (102), 
we arrive at the following result: 


Am2E2y2 Am? 
—_ Biers 
Ls aty2 i 5°. 

This constitutes the analytical description 
of Kepler’s third law: the squares of the 
periods of revolution of the planets are pro- 


portional to the cubes of the major azes of the 
planets’ orbits. 


1.19 Beam Deflection 


Let us consider a horizontal beam AB 
(Figure 32) of a constant cross section and 
made of a homogeneous material. The 
symmetry axis of the beam is indicated in 
Figure 32 by a dashed line. Suppose that 
forces acting on the beam in the vertical 
plane containing the symmetry axis bend 
the beam as shown in Figure 33. These forces 
may be the weight of the beam itself or 
an external force or the two forces acting 
simultaneously. Clearly, the symmetry axis 
will also bend due to the action of these 
forces. Usually a bent symmetry axis is 
called the elastic line. The problem of deter- 
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mining the shape of this line plays an im- 
portant role in elasticity theory. 

Note that there can be various types of 
beam depending on the way in which beams 
are fixed or supported. For example, Fig- 
ure 34 depicts a beam whose end A is rigidly 
fixed and end B is free. Such a beam is said 
to be a cantilever. Figure 35 depicts a beam 
lying freely on supports A and B. Another 
type of beam with supports is shown in 
Figure 36. There are also various ways in 
which loads can be applied to beams. For 
example, a uniformly distributed load is 
shown in Figure 34. Of course, the load can 
vary along the entire length of the beam 
or only a part of this length (Figure 35). 
Figure 36 illustrates the case of a concen- 
trated load. 

Let us consider a horizontal beam OA 
(Figure 37). Suppose that its symmetry 
axis (the dashed line in Figure 37) lies on 
the x axis, with the positive direction 
being to the right of point O, the origin. 
The positive direction of the y axis is 
downward from point O. External forces 
fas. T's; (and the weight of the beam 
if this is great) bend the symmetry axis, 
which becomes the elastic line (also depict- 
ed in Figure 38 by a curved dashed line). 
The deflection y of the elastic line from the 
x axis is known as the sag of the beam at 
point z. Thus, if we know the equation of 
the elastic line, we can always find the sag 


8—0770 
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of the beam. Below we show how this can 
be done in practical terms. 

Let us denote by M (x) the bending mo- 
ment in the cross section of the beam at 
coordinate xz. The bending moment is de- 
fined as the algebraic sum of the moments 
of the forces that act from one side of the 
beam at point zx. In calculating the mo- 
ments we assume that the forces acting on 
the beam upward result in negative mo- 
ments while those acting downward result 
in positive moments. 
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The strength-of-materials course proves 
that the bending moment at point z is re- 
lated to the curvature radius of the elastic 


line via the equation 


y” a 
EY Tye = M @) 


(103) 


where E is Young’s modulus, which de- 


8% 
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Fig. 38 


pends on the type of material of the beam, 
J is the moment of inertia of the cross sec- 
tion of the beam at point z about the hori- 
zontal straight line passing through the 
center of mass of the cross section. The 
product EJ is commonly known as the 
flexural rigidity, in what follows we assume 
this product constant. 

Now, if we suppose that the sag of the 
beam is small, which is usually the case in 
practice, the slope y’ of the elastic line is 
extremely small and, therefore, instead of 
Eq. (103) we can take the approximate 
equation 


EJy" = M (2). (104) 


To illustrate how Eq. (104) is used in 
practice, we consider the following prob- 
lem. A horizontal homogeneous steel beam 
of length / lies freely on two supports and 


= 
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Fig. 39 


sags under its own weight, which is p kgf 
per unit length. We wish to determine the 
equation of the elastic line and the maxi- 
mal sag. 

In Figure 39 the elastic line is depicted 
by the dashed curve. Since we are dealing 
with a two-support beam, each support acts 
on the beam with an upward reaction force 
equal to half the weight of the beam (or 
pl/2). The bending moment M (z) is the 
algebraic sum of the moments of these forces 
acting on one side from point Q (Fig- 
ure 39). Let us first consider the action of 
forces to the left of point Q. At a distance z 
from point Q a force equal to pl/2 acts on 
the beam upward and generates a negative 
moment. On the other hand, a force equal 
to px acts on the beam downward at a 
distance x/2 from point Q and generates a 
positive moment. Thus, the net bending 
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moment at point Q is given by the formula 


pl |. x px? plx 
M (2)= —-y- e+ ps (> )=4-— AP 


(105) 


If we consider the forces acting to the 
right of point Q, we find that a force equal 
to p (l — x) acts on the beam downward. at 
a distance (1 — x)/2 from point Q and gen- 
erates a positive moment, while a force 
equal to pl/2 acts on the beam upward at a 
distance | — x from point Q and generates 
a negative moment. The net bending mo- 
ment in this case is calculated by the for- 
mula 


[— l 
M (2) = p (l— x) —— — = (l—2) 
_ px® plz 
=“ 7S (106) 


Formulas (105) and (106) show that the 
bending moments prove to be equal. Now, 
knowing how to find a bending moment, we 
can easily write the basic equation (104), 
which in our case assumes the form 


Bi ipf a 22 222 (107) 


Since the beam does not bend at points O 
and A, we write the boundary conditions 
for Eq. (107) so as to find y: 


y=Oatz=—Oandy=Oatr=l. 
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Integration of Eq. (107) then yields 
Y = pep (x! — 2lz3 + Ba). (108) 


This constitutes the equation of the elastic 
jine. It is used to calculate the maximal 
sag. For instance, in this example, basing 
our reasoning on symmetry considerations 
(the same can be done via direct calcula- 
tions), we conclude that the maximal sag 
will occur at z= 1/2 and is equal to 
5pl*/384EJ, with EF = 21 x 10° kgf/cm? 
and J = 3 x 10% cm*. 


1.20 Transportation of Logs 


In transporting logs to saw mills, logging 
trucks move along forest roads some of the 
time. The width of the forest road is usual- 
ly such that only one truck can travel along 
the road. Sections of the road are made wid- 
er so that trucks can pass each other. Ignor- 
ing the question of how traffic should be 
arranged that loaded and empty traffic 
trucks meet only at such sections, let us 
establish how wide the turns in the road 
must be and what trajectory the driver 
must try to follow so that, say, thirty-me- 
ter logs can be transported. It is assumed 
that the truck is sufficiently maneuverable 
to cope with a limited section of the road. 

Usually a logging truck consists of a 
tractor unit and a trailer connected freely 
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Fig. 40 


to each other. The tractor unit has a front 
(driving) axle and two back axles above 
which a round platform carrying two posts 
and a rotating beam are fastened. This 
platform can rotate in the horizontal plane 
about a symmetrically positioned vertical 
axis. One end of each log is placed on this 
platform. The trailer has only two back 
axles, but also has a platform with two 
posts. The other ends of the logs are put on 
this platform. The trailer’s chassis consists 
of two metal cylinders one of which can 
partly move inside the other. The chassis 
connects the back platform with the axis 
that links the trailer with the tractor unit. 
Thus, the length of the chassis can change 
during motion, which enables the tractor 
unit and the trailer to move independently 
to a certain extent. Schematically the log- 
ging truck is depicted in Figure 40. The 
points A and B correspond to the axes of 
the platforms a distance h apart. By XY 
we denote a log for which AX = dh. Point 
C corresponds to a small axis that connects 
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the tractor unit with the trailer, with 
AC =ah. Usually a = 0.3, but in the 
simplest case of log transportation a = 0. 
Next, FF is the front axle of the tractor 
unit, and PP and QQ are the back axles of 
the tractor unit, while RR and SS are 
the trailer axles. All axles have the same 
length 2Z, so that for the sake of simplic- 
ity we assume that the width of the log- 
ging truck is 2L. As for the width of the 
load in its rear, we put it equal to 2W. 
In what follows we will need the concept 
of the sweep of the logging truck, which is 
commonly understood to be the maximum 
deflection of the rear part of the logging 
truck (for the sake of simplicity we assume 
this part to be point X) from the trajectory 
along which the logging truck moves. Let 
us suppose that the road has a width of 
2Bh and that usually a turn in the road is 
simply an arc of a circle of radius h/a 
centered at point O (Figure 41). For the 
sake of simplicity we assume that a logging 
truck enters a turn in the road in such a 
way that the tractor unit and the trailer 
are positioned along a single straight line, 
the driver operates the truck in such a way 
that point A, corresponding to the axis of 
the front platform, is exactly above the 
road’s center line. Point A in Figure 41 
is determined by the angle y that the truck 
AC makes with the initial direction. Here 
it is convenient to fix a coordinate system 
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Initial 
direction 


Fig. 41 


Oxy in such a way that the horizontal axis 
points in the initial direction and the ver- 
tical axis is perpendicular to it. In a general 
situation the load carried by the truck 
will make an angle 6 with the initial direc- 
tion. As for angle BAC in Figure 41, 
denoting it by wu we find that u = y — 0. 
Usually this angle is the logging truck’s 
angle of lag. The required halfwidth h of 
the road, which determines the sweep of the 
logging truck at a turn and is known as the 
halfwidth of the road at the outer curb of a 
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turn, is defined as the algebraic sum 
OX — OA + W, while the halfwidth of the 
road at the inner curb of a turn is defined 
as the algebraic sum OA + L — OP, where 
OP is the perpendicular dropped from 
point O onto AB. 

We stipulate that in its motion a logging 
truck’s wheel either experiences no lateral 
skidding at all or the skidding is small. 
This requirement, for one thing, means that 
the center line AC of the tractor unit con- 
stitutes a tangent to the arc of a circle at 
point A, so that OA is perpendicular to AC, 
and angle y is determined by the motion of 
point A along the arc of the circle. Note, 
further, that in building a road the curva- 
ture of a turn in the road is determined by 
an angle N° that corresponds to an arc 
length of a turn of approximately 30 m. 
In our notation, 

> 180° 30% 

N°= ae ae (109) 


where h is measured in meters. Thus, at 
h=9m and a = 0.1 we have N° w 19°, 
while at h = 12 m and a = 1.0 we have 
N° ~ 142°. For practical considerations we 
must consider only such a@’s that lie be- 
tween 0 and1,and the greater the value of 
a, the greater the maneuverability of the 
logging truck. 

The log length Ah will be greater than A, 
but again, reasoning on practical grounds, 
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it must not be greater than 3h. Thus, the 
value of A’ varies between 1 and 3. As for 
constant a, it is assumed that 0 < a < 0.5. 
Finally, we note that in each case the 
value of h is chosen differently, but it varies 
between 9 and 12 m. 

Since the wheels of the tractor unit do 
not skid laterally, the coordinates of point 
A in Figure 41 are 


hh. 
= — sin X, y = cos X. 
a a 


The coordinates of point B are 


X= sin X—hcos 8, 
Y= cos% +hsin0. (410) 


Since the trailer’s wheels do not skid either, 
point B moves in the direction BC and 


dy 


where wis the angle which BC makes with 
the initial direction. Next, employing the 
fact that y = wu + 6 and studying triangle 
ABC, we arrive at the following chain of 
equalities: 


sin(y—v) _ sin(0—) __ sinu 
h _— ah ~ bh ? (112) 


where 0< b< i, and wy, 0, and w are 
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functions of y. If we combine (411) with 
(110), we find that 


h. do 
(-— sin X +h Fy 008 8} cos > 
+-( = cos % +h sin @) sin p=0. 


Carrying out the necessary calculations, we 
arrive at 


00s (9 — wy). (113) 
If we exclude variable from this relation- 
ship for a fixed value of a by employing (112) 
and the fact that y =u 6, we arrive 
(since 6 = 0 at y = 0) at the differential 
equation 


du __y__ sin (114) 


dx a (1—a cos u) 


sin (X—) =a 


with the initial condition u (0) = 0, where 
the angle of lag plays the role of the sought 
function.* 

By substituting v for tan (u/2) in the 
differential equation (414), we can integrate 
the equation in closed form. However, the 
resulting wu vs. y dependence proves to be 
extremely complicated, which, of course, 
hinders an effective study. Nevertheless, 
Eq. (4414) can easily be integrated and 
studied numerically. To this end we can 


* See A.B. Tayler, “The sweep of a logging truck”, 
Math. Spectrum 7, No. 1: 19-26 (1974/75). 
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employ, for example, the Runge-Kutta sec- 
ond-order numerical method the essence of 
which is the following. 

Suppose that an integral curve u = @ (4) 
of a differential equation du/dy = f (y, u) 
passes through a point (yp, Wy). At equidis- 
tant points Yo, Xi Xo +> (Kita — Xi = 
Ay > 0) we select values uy, U4, Us, - 
such that u; ~ @; (y), where the successive 


values U,, Us, are specified by the 
formulas 

Wit, = Ui + (ky + k,)/2, 

with 


ky = i (X21 u;) AY, 
ke =f (yi + Ay, ui + hy) Ay. 


To solve Eq. (414) numerically with the 
initial condition u (0) = 0, we compile the 
following program using BASIC: 


140 REM Runge-Kutta second-order method 
20 DEF FNF (X, U) = 
4 — SIN(U)/(AL#(1 — A+COS(U))) 
30 CLS: PRINT “Solution of differential equa- 
tion” 
40 PRINT “Runge-Kutta second-order method” 
100 PRINT “Parameters:” 
1410 INPUT “Alpha=”, AL 
120 INPUT “a=”, A 
130 INPUT “Step of independent variable=”, DX 
140 PRINT “Initial values:” 
150 INPUT “of independent variable=”, X 


elise 
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160 INPUT “of function=”, U 

200 REM Next 20 values of function 
210 CLS: PRINT “X”, “U” 

220 FOR I = 1 TO 20 

230 PRINT X, U 

240 K1 = DX«*FNF (X,U) 

250 X = X + DX 

260 K2 = DX*FNF(X, U + K1) 

270 U = U+ (K1 4+ K2)/2 

280 NEXT I 

290 INPUT “Continue (Yes = 1, No = 0)”; I 
300 IFI(¢ ) 0 GOTO 210 

310 END 


Note that to compile the table of values 
of the solution for concrete values of a and 
a, we must know the limiting value C of 
this solution. The number C can be found 
from the condition that 


a(t—acosC) =sinC, 
which leads to the formula 


Cam sint (Se VE at taal) 


Here, if a=—1, we have = n/2 — 
2 tana, butifa < 1, thenC wa (1 — a). 
The limiting value C can be approximated 
by an exponential, namely 


C—unve-v%, 


where y = a (cos C — a)/sin? C. For small 
a’s the value of y is fairly great and is 
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approximately given by the relationship 


ne ae 
a (1 —a) 


Vr 


But if a = 1, we have 


{ 2 
ae 


For the step Ay in the variation of the 
independent variable y we must take a 
number that does not exceed C. This require- 
ment is especially important for small a’s 
and in the neighborhood of y = 0. Below 
we give the protocol for calculating the 
values of the function u = u (y). 


Solution of differential equation by 
Runge-Kutta second-order method 


Parameters: 

alpha = 1.0 

a= 0.3 

Step of independent variable = 0.5 
Initial values: 

of independent variable = 0 

of function = 0 


vs u 
0 0 

2 1437181 
4 .2299778 
6 .2824244 
8 .3146894 
1 3347111 
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2.4 

2.600001 
2.800001 
3.000001 
3.200001 
3.400001 
3.600001 
3.800001 


Continue (Yes = 1, No = 0)? 1 


xX 


4.000001 
4.200001 
4.4 

4.6 

4.8 

° 

0.2 
5.399999 
5.599999 
5.799999 
5.999999 
6.199999 
6.399998 
6.599998 
6.799998 
6.999998 


9—0770 


.3472088 
.3900403 
-30996 

.3630556 
.3650054 
.3662343 
.3670091 
.3674978 
.367806 

.3680005 
.3681 232 
-3682006 
.3682494 
.3682802 


u 


.3682997 
.368312 

.3683197 
.3683246 
.3683276 
.3683296 
.3683308 
-3683316 
.3683321 
.3683324 
.3683326 
.3683327 
.3683328 
.3683328 
.3683329 
.3683329 
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7.199998 .3683329 
7.399997 .3683329 
7.599997 .3683329 
7.799997 .3683329 
Continue (Yes = 1, No = 0)? 1 

x, u 
7.999997 .3683329 
8.199997 .3683329 
8.399997 .3683329 
8.599997 .3683329 
8.799996 .3683329 
8.999996 .3683329 
9.199996 .3683329 
9.399996 .3683329 
9.599996 .3683329 
9.799996 .3683329 
9.999995 .3683329 
10.2 .3683329 
10.4 .3683329 
10.6 .3683329 
10.8 .3683329 
10.99999 .3683329 
11.19999 .3683329 
11.39999 .3683329 
11.59999 .3683329 
11.79999 .3683329 
Continue (Yes = 1, No = 0)? 0 
OK 


The results of a numerical calculation at 
a =0Q.3 are presented graphically in Fig- 
ure 42. They show how an increase in ¥ 
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influences the angle of lag of the logging 
truck. For clarity of exposition the scales 
on the uw and y axes are chosen to be differ- 
ent. 

Now let us determine the sweep of the 
logging truck by using the concept of half- 
width of the road at the outer curb of a turn 
(this quantity was defined earlier as the al- 
gebraic sum OX — OA + W (Figure 41)). 
First we note that 


OX? = (— sin X— Ah cos 8)? 


+(+ cos X + An sin 6)? 
=: fy2 (+2224 sin u). =~ 


This shows that the sweep decreases as ¥ 
grows since the angle of lag, u, increases. 
Thus, the maximal halfwidth Bh of the 
road can be determined from the following 
formula for B: 


1 1 W 
B=V Mia ety 


In Figure 43 the solid curves reflect the 
relationship that exists between B — W/h 
and a for different values of A, while the 
dashed curves show what must be the 
value of 8 — L/h to guarantee the necessary 
“margin” at the inner curb of the turn. For 


Qe 
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C= 09878826 


C=0.5683329 
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C=0.0701512 
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every position of the logging truck on the 
road we must have 


B= = (4 — COS u) + 
<5 (1—cos €) +4 (115) 


Next, since the value of C decreases with a, 
the angle of lag u proves to be the greatest 
when a =0, which corresponds to the simpl- 
est case of logging. Then (115) yields 


1— V1i—a? 


L 1 
p— —-<— (1 —cos C) oe - 


Now let us dwell on the results that follow 
from the above line of reasoning. Here is a 
typical example that illustrates these re- 
sults. 

If a logging truck in which the distance 
between the front and back platforms is 
12 m follows a turn along an arc of a circle 
of radius 60 m, then @ = 0.2 and, in accor- 
dance with (109) N° is approximately 28°. 
If the width of the logging truck is 2.4 m 
and width of the load in the rear is 1.2 m, 
then for logs of approximately 24 m in 
length (reckoning from the axis of the 
front platform) we have A= 2, Wilh = 
0.05, and L/h = 0.1. Thus, from Figure 43 
it follows that for all values of a we have 
6B = 0.45 for the outer curb of the turn and 
B = 0.2 for the inner curb. Theoretically 
the necessary halfwidth of the road at the 
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outer curb of the turn is equal to 5.4 m 
and that at the inner curb to 2.4 m. [f the 
logging truck transports logs whose length 
is 14.4 m (reckoned from the axis of the 
front platform) and whose bunch width in 
the rear is 1.8 m, then in the case at hand 
X) = 1.2. The value of 6 then proves to be 
the same for the inner and outer curbs of 
the turn and equal to 0.22. Hence, as can 
easily be seen, the necessary halfwidth of 
the road at the outer curb of the turn is 
2.64 m, while the same quantity at the 
inner curb, as in the previous case, is equal 
to 2.4 m. This reasoning shows that the 
longer the logs transported the wider the 
road at a turn must be. For one, if we com- 
pare the two cases considered here, we see 
that an increase in the length of the logs by 
9.6 m requires widening the road at a turn 
by 2.76 m so that the driver can drive the 
truck along a curve whose length at the turn 
is approximately equal to the length of the 
road’s center line. Practice has shown that 
an inexperienced driver is not able to 
drive his truck along such a curve and needs 
a road whose total width at a turn is at 
least 10.8 m (if the load transported is 24 m 
long) for a truck width of 2.4 m. 

The theory developed above shows that 
the sweep of the logging truck is the great- 
est when the truck enters the turn, since 
in this case the angle of lag grows. This 
conclusion also holds true in the situation 
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when the truck enters one section of a zigzag 
turn after completing the previous one at 
the point of inflection. The results represent- 
ed graphically in Figure 43 correspond to 
the case where prior to entering a turn the 
tractor unit and the trailer are positioned 
on a single straight line. But if there exists 
a nonzero initial angle of lag C, caused by 
the zigzag nature of the turn, we must select 
the initial condition in the differential 
equation (114) in the form wu (0) = —C,. 
Then the required width of the road is 
determined from the following formula for B: 


sre | Ao: 4 W 
6 = V heap te = BIN Cy + ; 


We note that in the case of simple log 
transportation, that is,a=0, it is impossible 
to pass a turn with @ greater than unity. 
However, for fairly large values of a the 
values @ >> 1 become possible, they must 
obey the following relationship: a (14 — 
acos C) = sin C. Thus, the maximal value 
of a is (4 — a*)-!/*, with the practical 
extremal value of a being 1.25 at a = 0.5. 

We note in conclusion that for @ > 0.5 
considerable economy in the width of a road 
is achieved by increasing a (see Figure 43). 
If the load is such that A is not much greater 
than unity, the value of @ is chosen such 
that the required halfwidth of the road at 
the inner curb of any turn is always smaller 
than at the outer curb. 


Chapter 2 

Qualitative Methods 

of Studying Differential 
Models 


In solving the problems discussed in Chap- 
ter 1 we constructed differential models 
and then sought answers by integrating the 
resulting differential equations. However, 
as noted in the Preface, the overwhelming 
majority of differential equations are not 
integrable in closed (analytical) form. 
Hence, to study differential models of real 
phenomena and processes we need methods 
that will enable us to extract the necessary 
information from the properties of the 
differential equation proper. Below with 
concrete examples we show how in solving 
practical problems one can use the simplest 
approaches and methods of the qualitative 
theory of ordinary differential equations. 


2.4 Curves of Constant Direction 
of Magnetic Needle 


Let us see how in qualitative integration, 
that is, the process of establishing the gen- 
eral nature of solutions to ordinary differen- 
tial equations, one can use a general proper- 
ty of such equations whose analogue, for 
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example, is the property of the magnetic 
field existing at Earth’s surface. The reader 
will recall that curves at Earth’s surface 
can be specified along which the direction 
of a magnetic needle is constant. 

Thus, let us consider the first-order or- 
dinary differential equation 


SU f(a, y), (116) 


where the function f (z, y) is assumed single- 
valued and continuous over the set of 
variables x and y within a certain domain D 
of the (z, y)-plane. To each point M (z, y) 
belonging to the domain D of function 
f (x, y) the differential equation assigns 
a value of dy/dz, the slope K of the tangent 
to the integral curve at point M (z, y). 
Bearing this in mind, we say that at each 
point M (za, y) of D the differential equation 
(416) defines a direction or a line element. 
The collection of all line elements in D is 
said to be the field of directions or the 
line element field. Graphically a linear ele- 
ment is depicted by a segment for which 
point M (z, y) is an interior point and which 
makes an angle 0 with the positive direction 
of the x axis such that K =tan 0= 
f (xz, y). From this it follows that geometri- 
cally the differential equation (116) expresses 
the fact that the direction of a tangent at 
every point of an integral curve coincides with 
the direction of the field at this point, 
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To construct the field of directions it has 
proved expedient to use the concept of iso- 
clinals (derived from the Greek words for 
“equal” and “sloping”), that is, the set of 
points in the (2, y)-plane at which the 
direction of the field specified by the differ- 
ential equation (116) is the same. 

The isoclinals of the magnetic field at 
Earth’s surface are curves at each point of 
which a magnetic needle points in the same 
direction. As for the differential equation 
(116), its isoclinals are given by the equa- 
tion 
f (x, y) = %, 


where v is a varying real parameter. 

Knowing the isoclinals, we can approx- 
imately establish the behavior of the 
integral curves of a given differential equa- 
tion. Let us consider, for example, the 
differential equation 


dy __ 2 2 
Pe + Yy*, 


which cannot be integrated in closed (ana- 
lytical) form. The form of this equation 
suggests that the family of isoclinals is given 
by the equation 


ety=v, v>), 


that is, the isoclinals are concentric circles 


of radius V v centered at the origin and ly- 
ing in the (z, y)-plane, At each point of 


= 
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Fig. 44 


such an isoclinal the slope of the tangent 
to the integral curve that passes through 
this point is equal to the squared radius of 
the corresponding circle. This information 
alone is sufficient to convey an idea of the 
behavior of the integral curves of the 
given differential equation (Figure 44). 

We arrived at the final result quickly 
because the example was fairly simple. 
However, even with more complicated equa- 
tions knowing the isoclinals may prove to 
be expedient in solving a specific problem. 

Let us consider a geometric method of 
integration of differential equations of the 
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type (116). The method is based on using 
the geometric properties of the curves 
given by ha equations 


f(e, y)= (0) 
of ee y) oie Ane 7 _ Y) f(z, y)=0. (L) 


Equation (J - a is, the “zero” isoclinal 
equation, specifies curves at whose points 
dy/dz = 0. This means that the points of 
these curves may prove to be points of maxi- 
ma or minima for the integral curves of 
the initial differential equation. This ex- 
plains why out of the entire set of isoclinals 
we isolate the “zero” isoclinal. 

For greater precision in constructing 
integral curves it is common to find the 
set of inflection points of these curves (pro- 
vided that such points exist). As is known, 
points of inflection should be sought among 
the points at which y” vanishes. Employing 
Eq. (416), we find that 


n __ Of(z, Of (z, 
y’ = ARSON LE STE Wy! 


__ Of (z, y) Of (z, y) 
=A | SEY) Fe, y). 


The curves specified by Eq. (L) are the pos- 
sible point-of-inflection curves.* Note, for 


* It is assumed here that integral curves that 
fill a certain domain possess the property that only 
one integral curve passes through each point of the 
domain, 
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one, that a point of inflection of an integral 
curve is a point at which the integral curve 
touches an isoclinal. 

Thezcurves consisting of extremum points 
(maxima and minima points) and points of 
inflection of integral curves break down the 
domain of f into such subdomains S,, 
So, +«, Sm in which the first and second 
derivatives of the solution to the differ- 
ential equation have definite signs. In 
each specific case these subdomains should 
be found. This enables giving a rough 
picture of the behavior of integral curves. 

As an example let us consider the differ- 
ential equation y’ = x -+ y. The equation 
of curve (J,) in this case has the form 
x+y=0, or y=-—z. A direct check 
verifies that curve (J,) is not an integral 
curve. As for curve (Z), whose equation in 
this case is y+2-+1=0, we find that 
it is an integral curve and, hence, is not 
a point-of-inflection curve. 

The straight lines (J,) and (ZL) break 
down the (z, y)-plane into three subdomains 
(Figure 45): S, (y’ >0, y” >0), to the 
right of the straight line (J,), S. (y’ <0, 
y” > 0), between the straight lines (J,) and 
(L), and S;,(y’< 0, y” <0), to the left of the 
straight line (Z). The points of minima of the 
integral curves lie on the straight line (/,). 
To the right of (7,) the integral curves point 
upward, while to the left they point down- 
ward (left to right in Figure 45). There are 
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no points of inflection. To the right of the 
straight line (Z) the curves are convex down- 
ward and-to the left convex upward. The 
behavior of the integral curves on the whole 
is shown in Figure 45. 

Note that in the given case the integral 
curve (L) is a kind of “dividing” line, since 
it separates one family of integral curves 
from another. Such a curve is commonly 
known as a separatrix. 
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2.2 Why Must an Engineer Know 
Existence and Uniqueness 
Theorems? 


When speaking of isoclinals and point-of- 
inflection curves in Section 2.1, we tacitly 
assumed that the differential equation in 
question had a solution. The problem of 
when a solution exists and of when it is 
unique is solved by the so-called existence 
and uniqueness theorems. These theorems 
are important for both theory and practice. 

Existence and uniqueness theorems are 
highly important because they guarantee 
the legitimacy of using the qualitative meth- 
ods of the theory of differential equations 
to solve problems that emerge in science 
and engineering. They serve as a basis for 
creating new methods and theories. Often 
their proof is constructive, that is, the meth- 
ods by which the theorems are proved sug- 
gest methods of finding approximate solu- 
tions with any degree of accuracy. Thus, 
existence and uniqueness theorems lie at 
the base of not only the above-noted quali- 
tative theory of differential equations but 
also the methods of numerical integration. 

Many methods of numerical solution of 
differential equations have been developed, 
and although they have the common draw- 
back that each provides only a concrete 
solution, which narrows their practical 
potential, they are widely used. It must be 
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noted, however, that before numerically 
integrating a differential equation one must 
always turn to existence and unique- 
ness theorems. This is essential to avoid 
misunderstandings and incorrect conclu- 
sions. 

To illustrate what has been said, let us 
take two simple examples,* but first let us 
formulate one variant of existence and 
uniqueness theorems. 

Existence theorem If in Eq. (416) func- 
tion f is defined and continuous on a bounded 
domain D in the (x, y)-plane, then for every 
point (29, Yo) € D there exists a solution y (zx) 
to the initial-value problem ** 


d 
se =f le, y)s YY (20s) = Yor (117) 


that is defined on a certain interval con- 
taining point Zz. 

Existence and uniqueness theorem If 
in Eq. (A416) function f is defined and contin- 
uous on a bounded domain D in the (z, y)- 


* See C.E. Roberts, Jr., “Why teach existence and 
uniqueness theorems in the first course of ordinary 
differential equations?”, Int. J. Math. Educ. Sci. 
Technol. 7, No. 1: 41-44 (1976). 

** If we wish to find a solution ofa differential 
equation satisfying a certain initial condition (in 
our case the initial condition is y (x9) = yy), such 
a problem is said to be an initial-value problem. 
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plane and satisfies in D the Lipschitz condi- 
tion in variable y, that is, 


lf @, v2) —f, wIl<Llye—y I, 


with L a positive constant, then for every 
point (to, Yo) € D there exists a unique solu- 
tion y (x) to the initial-value problem (117) 
defined on a certain interval containing 
point Lo. 

Extension theorem IJf the hypotheses of 
the existence theorem or the existence and 
uniqueness theorem are satisfied, then every 
solution to Eq. (A416) with the initial data 
(29, Yo) € D can be extended to a point that 
lies as close to the boundary of D as desired. 
In the first case the extension is not necessarily 
unique while in the second it is. 

Let us consider the following problem. 
Using the numerical Euler integration 
method with the iteration scheme y;4, = 
y; + hf (z;, yi) and step h = 0.4, solve 
the initial-value problem 


y’ = —aly, y(—1) = 0.24 (418) 


on the interval [— 4, 3]. 

Note that the problem involving the 
equation y’ = —x/y emerges, for example, 
from the problem considered on p. 162 con- 
cerned with a conservative system consist- 
ing of an object oscillating horizontally in 


10—0770 
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a vacuum under forces exerted by linear 
springs. 

To numerically integrate the initial-value 
problem (117) on the [—1, 2.8] interval 
we compile a program for building the 
graph of the solution: 


40 REM Numerical integration of differential 
equation 

20 DEF FNF(X, Y) = 

30 GOSUB 1110: REM Coordinate axes 

40 REM Next value of function 

1000 LINE —(FNX(X), FNY(Y)), 2 

1010 IF X < 2.9 GOTO 100 

1020 LOCATE 23, 14 

1030 END 

1100 REM Construction of coordinate axes 

1110 SCREEN 1, 1, 0: KEY OFF: CLS 

1120 DEF FNX(X) = 88 + 80*X 

1130 DEF FNY(Y) = 96 — 80*Y 

1140 REM Legends on coordinate axes 

1150 LOCATE 4, 11, 0 PRINT “Y” 

4160 LOCATE 3, 14: PRINT “4” 

1170 LOCATE 13, 39: PRINT “X” 

11480 LOCATE 23, 10: PRINT “—1” 

11490 FOR I = —1 TO 2: LOCATE 143, 10«I + 
41: PRINT USING “# +”: I;:; NEXT I 

1200 REM Oy 

1210 DRAW “BM88, ONM — 2, +8NM -+ 2, 
+8D16” 

1220 FORI = 1T011: DRAW “NR2D16": NEXT 

1230 REM Ox 
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1240 DRAW “BM319, 96NM — 7, —2NM — 7 
+2123” 

1250 FORI = 1 TO 19: DRAW “NU2L16”: NEXT 

4260 REM Input of initial data 

1270 LOCATE 25, 1: PRINT STRING$(40, —_”), 

4280 LOCATE 25, 1: INPUT, “xO =”, X: INPUT; 
“dx =”, DX: INPUT; “yO =”, Y 

1290 PSET (FNX(X), FNY(Y)), 2 

1300 RETURN 


Here in line 20 the right-hand side of the 
differential equation is specified, and lines 
50 to 990 must hold the program of the 
numerical integration of the differential 
equation. 

In our case (the initial-value problem 
(1418)) the beginning of the program has 
the following form: 


10 REM Euler’s method 

20 DEF FNF(X, Y) = —X/Y 

30 GOSUB 1100: REM Coordinate axes 
40 REM Next value of function 

100 Y = Y + DX*FNF(X, Y) 

1440 X = X+ DX 


The results are presented graphically in 
Figure 46. 

Let us now turn to the existence theorem. 
For the initial-value problem (418), the 
function f(z, y) = —z/y is defined and 
continuous in the entire (z, y)-plane exclud- 
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ing the zx axis. Thus, in accordance with 
the existence theorem, the initial-value 
problem (118) has a solution y (z) defined on 
a certain interval containing point z, = 
—1. This solution, according to the exten- 
sion theorem, can be extended to a value 
of y (x) close to the value y (7) = 0. As 
a result of numerical integration we have 
arrived at a solution of (118) defined on an 
interval (a, b), with a< —1 and 1.3< 
b<1.4. However, allowing for the concrete 
form of the differential equation, we can 
specify the true interval in which the 
solution to the initial-value problem (418) 
exists. Indeed, since in the initial equation 
the variables can be separated, we have 


7] x 


ndy= — E dE. 
0.21 -i1 


Integrating, we get y = V 1.0441 — 2°. 
Hence, a solution to the initial-value prob- 
lem (448) exists only for | z |< VY 1.0441 
~ 1.0218. 

Thus, by resorting to the existence theo- 
rem (and to the extension theorem) we were 
able to “cut off” the segment on which there 
is certain to be no solution of the initial- 
value problem. If we employ only numerical 
integration, we arrive at an _ erroneous 
result. The fact is that as the solution 
y = y (x) approaches the z axis, the angle 
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of slope of the curve tends to 90°. Therefore, 
in the time that the independent varia- 
ble z changes by 0.1, the value of y is able to 
“jump over” the z axis, and we find ourselves 
on an integral curve that differs from the 
original. This happens because the Euler 
method allows for the angle of slope only 
at the running point. 

The following example is even more 
instructive. We wish to solve the initial- 
value problem 


y’ = 327y, y(—1) =—1 (119) 


on the segment [— 1, 1]. The approach here 
consists in first employing the Euler method 
and then an improved Euler method with 
a step h = 0.1 and an iteration scheme 
Vier = Yi thf Gita Yi+1/2)s with 
Yissso = Yi + hf (ai, yi)/2. 

We solve the initial-value problem (119) 
using the above program, whose beginning 
in the case at hand has the following form: 


10 REM Euler’s method 

20 DEF FNF(X, Y) =3*«X*SGN(Y)*ABS(Y) ~ 
(1/3) 

30 GOSUB 1100: REM Coordinate axes 

40 REM Next value of function 

100 Y = Y + DX*FNF(X, Y) 

410 X = X-+ DX 


The results are presented graphically in 
Figure 47. 


= 
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Fig. 47 


As for the improved Euler method, the 
beginning of the program has the following 
form: 


10 REM Improved Euler’s method 

20 DEF FNF(X, Y) = 3*X#SGN(Y)*ABS(Y) ~ 
(1/3) 

30 GOSUB 1100: REM Coordinate axes 

40 REM Next value of function 

50 D2 = DX/2 
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Fig. 48 


100 Yi = Y + D2«FNF(X, Y) 
140 Y = Y + DX«*FNF(X + D2, Y1) 
1200 X =X + Dx 


The results are presented graphically in 
Figure 48, and the diagram differs from 
that depicted in Figure 47. 

To look into the reason for such striking 
discrepancy in the results, we integrate the 
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initial-value problem. Separating the vari- 
ables, we get 


y x 
ve e = 3 ’ 
\n u3dy J Ee 


or, finally, y = +2°. This result already 
suggests that the solution via the Euler 
method gives the function y, (x) = z® while 
the solution via the improved Euler method 
gives 


| se oif «<0, 
Yo (x) = | —x? if xr>O0. 


Both y, and y, are solutions to the initial- 
value problem (119), which means that the 
solution of the initial-value problem con- 
sidered on the segment [—1, 1] is not 
unique. 

Let us now turn to the existence and 
uniqueness theorem in connection with this 
problem. First, we note that since the func- 
tion f (xz, y) = 3zVy is continuous in the 
entire (x, y)-plane, the existence theorem 
implies that the initial-value problem (119) 
has a solution defined on a segment con- 
taining point z, = —1, and, according to 
the extension theorem, this solution can 
be extended to any segment. Further, since 
Of (x, y)/dy = zy~-?/?, the function f (z, y)= 
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3xi/ y satisfies the Lipschitz condition in 
variable y in any domain not containing 
the x axis. If, however, a domain does 
contain points belonging to the z axis, it 
is easy to show that the function does 
not satisfy the Lipschitz condition. Hence, 
from the existence and uniqueness theorem 
(and the extension theorem) it follows that 
in this case the solution to the initial-value 
problem can be extended in a unique man- 
ner at least to the x axis. But since the 
straight line y = 0 constitutes a singular 
integral curve of the differential equation 


y’ = 3z%/y, we already know that as soon 
as y vanishes, there is no way in which we 
can extend the solution to the initial-value 
problem (419) beyond point O (0, 0) in 
a unique manner. 

Thus, by resorting to the existence and 
uniqueness theorem (and the extension 
theorem) we were able to understand the 
results of numerical integration, that is, if 
we are speaking of the uniqueness of the 
solution of the initial-value problem (119) 
on the [—1, 1] segment, the solution exists 
and is defined only on the segment [—1, 0]. 
Generally, however, there can be _ several 
such solutions. 
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2.3 A Dynamical Interpretation 
of Second-Order Differential 
Equations 


Let us consider the nonlinear differential 
equation 
d2z { dz 


rT eae ar a (120) 


whose particular case is the second-order dif- 
ferential equation obtained on p. 72 when 
we considered pendulum clocks. We take 
a simple dynamical system consisting of 
a particle of unit mass that moves along the 
x axis (Figure 49) and on which a force 
f (z, da/dt) acts. Then the differential 
equation (120) is the equation of motion 
of the particle. The values of x and dz/dt 
at each moment in time characterize the 
state of the system and correspond to a 
point in the (z, dz/dt)-plane (Figure 50), 
which is known as the plane of states or 
the phase (z, dz/dt)-plane. The phase plane 
depicts the set of all possible states of the 
dynamical system considered. Each new state 
of the system corresponds to a new point in 
the phase plane. Thus, the changes in the 
state of the system can be represented by 
the motion of a certain point in the phase 
plane. This point is called a representative 
point, the trajectory of the representative 
point is known as the phase trajectory, 
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Fig. 49 


0 x 
Fig. 50 


and the rate of motion of this point as the 
phase velocity. 

If we introduce the variable y = dz/dé, 
Eq. (420) can be reduced to a system of 
two differential equations: 


dz a 
ape — ap a, y). (121) 


If we take ¢ as a parameter, then the solu- 
tion to system (421) consists of two func- 
tions, x (t) and y (t), that in the phase 
(x, y)-plane define a curve (a phase trajec- 
tory). 
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It can be shown that system (121) and 
even a more general system 
dz d 
a= X(t, Y)s TP =Y (, y); (122) 
where the functions X and Y and their 
partial derivatives are continuous in a do- 
main D, posseses the property that if x (t) 
and y (t) constitute a’solution to the differ- 
ential system, we can write 


r=x(t+CcC)y=y(t+0), (423) 


where C is an arbitrary real constant, and 
(123) also constitute a solution to the same 
differential system. All solutions (123) 
with different values of C correspond to a 
single phase trajectory in the phase (z, y)- 
plane. Further, if two phase trajectories 
have at least one common point, they 
coincide. Here the increase or decrease in 
parameter ¢ corresponds to a certain direc- 
tion of motion of the representative point 
along the trajectory. In other words, a 
phase trajectory is a directed, or oriented, 
curve. When we are interested in the 
direction of the curve, we depict the direc- 
tion of the representative point along the 
trajectory by placing a small arrow on the 
curve. 

Systems of the (122) type belong to the 
class of autonomous systems of differential 
equations, that is, systems of ordinary 
differential equations whose right-hand 
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sides do not explicitly depend on time ¢. 
But if at least in one of the equations of 
the system the right-hand side depends 
explicitly on time ¢, then such a system is 
said to be nonautonomous. 

In connection with this classification of 
differential equations the following remark 
is in order. If a solution z (¢) to Eq. (420) 
is a nonconstant periodic solution, then 
the phase trajectory of the representative 
point in the phase (z, y)-plane is a simple 
closed curve, that is, a closed curve without 
self-intersections. The converse is also true. 

If differential systems of the (122) type 
are specified in the entire (x, y)-plane, then, 
generally speaking, phase trajectories will 
completely cover the phase plane without 
intersections. And if it so happens that 


X (Xo, Yo) = Y (2, Yo) = O 


at a point M, (xo, yo), the trajectory degen- 
erates into a point. Such points are called 
singular. In what follows we consider primar- 
ily only isolated singular points. A singu- 
lar point M, (zo, Yo) is said to be isolated 
if there exists a neighborhood of this point 
which contains no other singular points 
except My (Xo, Yo). 

From the viewpoint of a physical inter- 
pretation of Eq. (120), the point M, (zo, 0) 
is a singular point. At this point y =0 and 
f (zo, 0) = O. Thus, in this case the isolated 
singular point corresponds to the state of 
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a particle of unit mass in which both the 
speed dz/dt and the acceleration dy/dt = 
d?x/dt? of the particle are simultaneously ze- 
ro, which simply means that the particle 
is in the state of rest or in equilibrium. 
In view of this, singular points are also 
called points of rest or points of equilibrium. 

The equilibrium states of a physical sys- 
tem constitute very special states of the 
system. Hence, a study of the types of 
singular points occupies an important place 
in the theory of differential equations. 

The first to consider in detail the clas- 
sification of singular points of differential 
systems of the (122) type was the distin- 
guished Russian scientist Nikolai BE. Zhu- 
kovsky (1847-1921). In his master’s the- 
sis “The kinematics of a liquid body”, 
presented in 1876, this problem emerged in 
connection with the theory of velocities 
and accelerations of fluids. The modern 
names of various types of singular points were 
suggested by the great French mathemati- 
cian Jules H. Poincaré (1854-1912). 

Now let us try to answer the question of 
the physical meaning that can be attached 
to phase trajectories and singular points of 
differential systems of the (422) type. For 
the sake of clarity we introduce a_two- 
dimensional vector field (Figure 51) defined 
by the function 


V (z, y) = X @, yi t+ Y (, y)j, 
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Fig. 54 


where i and j are the unit vectors directed 
along the x and y axes, respectively, in 
a Cartesian system of coordinates. At every 
point P (x, y) the field has two components, 
the horizontal X (z, y) and the vertical 
Y (2, y). Since da/dt = X (2, y) and dy/dt= 
Y (a, y), the vector associated with each 
nonsingular point P (x, y) is tangent at 
this point to a phase trajectory. 

If variable ¢ is interpreted as time, vec- 
tor V can be thought of as the vector of 
velocity of a representative point moving 
along a trajectory. Thus, we can assume 
that the entire phase plane is filled with 
representative points and that each phase 
trajectory constitutes the trace of a moving 
representative point. As a result we arrive 
at an analogy with the two-dimensional 
motion of an incompressible fluid. Here, 
since system (122) is autonomous, vector V 
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at each fixed point P (z, y) is time-inde- 
pendent and, therefore, the motion of the 
fluid is steady-state. The phase trajectories 
in this case are simply the trajectories of 
the moving particles of fluid and the singu- 
lar points O, O’, and O” (see Figure 51) are 
those where the fluid is at rest. 

The most characteristic features of the 
fluid motion shown in Figure 51 are (1) 
the presence of singular points, (2) the 
different patterns of phase trajectories near 
singular points, (3) the stability or instabil- 
ity at singular points (i.e. two possibili- 
ties may realize themselves: the particles 
that are in the vicinity of singular points 
remain there with the passage of time or 
they leave the vicinity for other parts of 
the plane), and (4) the presence of closed 
trajectories, which in the given case corre- 
spond .to periodic motion. 

These features constitute the main part 
of the phase portrait, or the complete 
qualitative behavior pattern of the phase 
trajectories of a general-type system (122). 
Since, as noted earlier, differential equa- 
tionscannot generally be solved analytically, 
the aim of the qualitative theory of ordinary 
differential equations of the (122) type is to 
build a phase portrait as complete as possi- 
ble directly from the functions X (z, y) 
and Y (z, y). 
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2.4 Conservative Systems 
in Mechanics 


Practice gives us ample examples of the 
fact that any real dynamical system 
dissipates energy. The dissipation usually 
occurs as a result of some form of friction. 
But in some cases it is so slow that it can be 
neglected if the system is studied over a fair- 
ly small time interval. The law of energy 
conservation, namely, that the sum of kinet- 
ic and potential energies remains constant, 
can be assumed to hold true for such sys- 
tems. Systems of this kind are called con- 
servative. For example, rotating Karth may 
be seen as a conservative system if we take 
a time interval of several centuries. But if 
we study Earth’s motion over several mil- 
lion years, we must allow for energy dissipa- 
tion related to tidal flows of water in seas 
and oceans. 

A simple example of a conservative sys- 
tem is one consisting of an object moving 
horizontally in a vacuum under forces ex- 
erted by two springs (Figure 52). If x is the 
displacement of the object (mass m) from 
the state of equilibrium and the force with 
which the two springs act on the object (the 
restoring force) is proportional to z, the 
equation of motion has the form 


d@x 


= 
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Fig. 52 


Springs of this type are known as linear, 
since the restoring force exerted by them is 
a linear function of z. 

If an object of mass m moves in a medi- 
um that exerts a drag on it (the damping 
force) proportional to the object’s velocity, 
the equation of motion for such a nonconser- 
vative — is 


m Se eo ~+4hkr=0, c> 0. (124) 


Here we are dealing with linear damping, 
since the damping force is a linear function 
of velocity dz/dt. 

If f and g are such arbitrary functions that 
f (0) = 0 and g (0) = O, the more general 
— 


a te (ar) t+f@=0 (125) 


can be interpreted as the equation of 
motion of an object of mass m under 
a restoring force —f (rz) and a damping 
force —g(dz/dt). Generally, these forces 
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are nonlinear; hence Eq. (125) can be con- 
sidered the basic equation of nonlinear me- 
chanics. 

Let us briefly examine the special case of 
a nonlinear conservative system described 
by the equation 

2x 

me +H (2)=0, (126) 
where the damping force is zero and, hence, 
energy is not dissipated. From Eq. (426) we 
can pass to the autonomous system 


dx dys f (2) 
Geet ap aa (127) 


If we now exclude time t from Eqs. (127), 
we arrive at a difierential equation for the 
trajectory of the system in the phase plane: 


ov a — Le) (1.28) 


dz my ° 
This equation can be written as 
my dy = —f (zx) dz. (129) 


Then, assuming that 7 = x, at ¢ = ¢, and 
y = Yo, We can integrate Eq. (129) from ¢, 
to ¢t. The result is 


x 
1 m 1 
+ myp—— my= — \ f( ak, 
Xo 
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which may be rewritten as 
x Xo 


s my? + | 7 © de= myt | f(x) de. 
‘ (130) 


Note that my?/2 = m (dz/dt)?/2 is the ki- 
netic energy of a dynamical system and 


x 


V (a) = | f@ ak (134) 


0 


is the system’s potential energy. Thus, 
Eq. (130) expresses the law of energy conser- 
vation: 


+ my? +V (2) =£, (132) 


where EF = my?/2 -++ V (a) is the total ener- 
gy of the system. Clearly, Eq. (132) is the 
equation of the phase trajectories of sys- 
tem (127), since it is obtained by integrating 
Eq. (128). Thus, different values of E corre- 
spond to different curves of constant energy 
in the phase plane. The singular points of 
system (127) are the points M, (zx,, 0), 
where x, are the roots of the equation f (z) = 
Q. As noted earlier, the singular points are 
points of equilibrium of the dynamical sys- 
tem described by Eq. (126). Equation (128) 
implies that the phase trajectories of the 
system intersect the z axis at right angles, 
while the straight lines z = zx, are parallel 
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to the x axis. In addition, Eq. (132) shows 
that the phase trajectories are symmetric 
with respect to the z axis. 

In this case, if we write Eq. (432) in the 
form 


y= +/ 21 -Via, (133) 


we can easily plot the phase trajectories. 
Indeed, let us introduce the (z, z)-plane, 
the plane of energy balance (Figure 53), 
with the z axis lying on the same vertical 
line as they axis of the phase plane. We then 
plot the graph of the function z = V (z) and 
several straight lines z = F in the (z, 2z)- 
plane (one such straight line is depicted in 
Figure 53). We mark a value of E — V (z) 
on the graph. Then for a definite x we multi- 
ply & — V (a) by 2/m and allow for formu- 
Ja (133). This enables us to mark the respec- 
tive values of y in the phase plane. Note 
that since dz/dt = y, the positive direction 
along any trajectory is determined by the 
motion of the representative point from left 
to right above the z axis and from right to 
left below the zx axis. 

The above reasoning is fairly general and 
makes it possible to investigate the equa- 
tion of motion of a pendulum ina medium 
without drag, which has the form (see p. 72) 


ar tksina=0, (134) 


= 
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Fig. 53 


with k a positive constant. 

Since Eq. (134) constitutes a particular 
case of Eq. (126), it can be interpreted as 
an equation describing the frictionless mo- 
tion in a straight line of an object of unit 
mass under a restoring force equal to 
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—k sin x. In this case the autonomous 
system corresponding to Eq. (134) is 


d d : 
=y, += —ksin x. (135) 


The singular points are (0, 0), (+a, 0), 
(+ 2n, 0), .., and the differential equa- 
tion of the phase trajectories of system (135) 
assumes the form 


dy __—s—s—sK sing 
dx y 


Separating the variables and integrating, 
we arrive at an equation for the phase tra- 
jectories, 


5 yk (1—cos z) =E. 


This equation is a particular case of Eq. (132) 
with m= 1, where the potential energy 
determined by (131) is specified by the 
relationship 


V (x) = ( f(t) dE =k(4—cosz). 


0 


In the (x, z)-plane we plot the function 
z == V (x) as well as several straight lines 
z= F (in Figure 54 only one such line, 
z = EF = 2k, is shown). After determining 
a value of E — V (z) we can draw a sketch 
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Fig. 54 


of the trajectories in the phase plane by 
employing the relationship 


y=+YV 2(E—V (x)}. 


170 Differential Equations in Applications 


The resulting phase portrait shows (see 
Figure 04) that if the energy £ varies from 0 
to 2k, the corresponding phase trajectories 
prove to be closed and Eq. (134) acquires 
periodic solutions. On the other hand, if 
E > 2k, the respective phase trajectories 
are not closed and Eq. (134) has no periodic 
solutions. Finally, the value EF = 2k corre- 
sponds to a phase trajectory in the phase 
plane that separates two types of motion, 
that is, is a separatriz. The wavy lines 
lying outside the separatrices correspond to 
rotations of a pendulum, while the closed 
trajectories lying inside the regions bound- 
ed by separatrices correspond to oscilla- 
tions of the pendulum. 

Figure 54 shows that in the vicinity of 
the singular points (+2nm, 0), m = 0, 1, 
2, .., the behavior of the phase trajecto- 
ries differs from that of the phase trajecto- 
ries in the vicinity of the singular points 
(nn, 0), wheren =1, 2,. .. 

There are different types of singular points. 
With some we will get acquainted shortly. 
As for the above example, the singular 


points (+ 2nm, 0), m =0, 1, 2, ..., be 
long to the vortez-point type, while the 
singular points (+n, 0), n=1, 2, .., 


belong to the saddle-point type. A singular 
point of an autonomous differential system 
of the (122) type is said to be a vortex 
point if there exists a neighborhood of this 
point completely filled with nonintersecting 
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phase trajectories surrounding the point. 
A saddle point is a singular point adjoined 
by a finite number of phase trajectories 
(“whiskers”) separating a neighborhood of 
the singular point into regions where the 
trajectories behave like a family of hyper- 
bolas defined by the equation zy = const. 

Now let us establish the effect of linear 
friction on the behavior of the. phase tra- 
jectories of a conservative system. The 
equation is 


3 
oe +e + ksine=0, c>0. 


If friction is low, that is, the pendulum is 
able to oscillate about the position of 
equilibrium, it can be shown that the phase 
trajectories are such as shown in Figure 50. 
But if friction is so high that oscillations 
become impossible, the pattern of phase 
trajectories resembles the one depicted in 
Figure 956. 

If we now compare the phase portrait of 
a conservative system with the last two por- 
traits of nonconservative systems, we see 
that saddle points have not changed (we 
consider only smal! neighborhoods of singu- 
lar points), while in the neighborhood of 
the singular points (+2am, 0), m = 0, 1, 
2, .., the closed phase trajectories have 
transformed into spirals (for low friction) 
or into trajectories that “enter” the singular 
points in certain directions (for high fric- 


R 
oy 
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tion). In the first case (spirals) we have 
singular points of the focal-point type and 
in the second, of the nodal-point type. 

A singular point of a two-dimensional 
autonomous differential system of the gen- 
eral type (122) (if such a point exists) 
is said to be a focal point if there exists 
a neighborhood of this point. that is com- 
pletely filled with nonintersecting phase 
trajectories resembling spirals that “wind” 
onto the singular point either as t—~ -+ co 
or as t-> —oo. A nodal point is a singular 
point in whose neighborhood each phase 
trajectory behaves like a branch of a parabo- 
la or a half-line adjoining the point along 
a certain direction. 

Note that if a conservative system has 
a periodic solution, the solution cannot be 
isolated. More than that, if [ is a closed 
phase trajectory corresponding to a periodic 
solution of the conservative system, there 
exists a certain neighborhood of I that is 
completely filled with closed phase tra- 
jectories. 

Note, in addition, that the above defi- 
nitions of types of singular points have 
a purely qualitative, descriptive nature. 
As for the analytical features of these types, 
there are no criteria, unfortunately, in 
the general case of systems of the (122) 
types, but for some classes of differential 
equations such criteria can be formulated. 
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dz d 

Fro = Ut + Oy, Fe = Ayt + bey, 


where a@,, 6,, a@,, and by, are real constants. 
If the coefficient matrix of this sytem is 
nonsingular, that is, the determinant 


a, Ob, 


#0, 


Ay Oe 


the origin O (0, 0) of the phase plane is 
the only singular point of the differential 
system. 

Assuming the last inequality valid, we de- 
note the eigenvalues of the coefficient ma- 
trix by A, and A,. It can then be demonstrat- 
ed that 

(1) if A, and A, are real and of the same 
sign, the singular point is a nodal point, 

(2) if A, and A, are real and of opposite 
sign, the singular point is a saddle point, 

(3) if A, and A, are not real and are not 
pure imaginary, the singular point is a fo- 
cal point, and 

(A) if A, and A, are pure imaginary, the 
singular point is a vortex point. 

Note that the first three types of singular 
points belong to the so-called coarse singu- 
lar points, that is, singular points whose 
nature is not affected by small perturhba- 
tions of the right-hand sides of the initial 
differential system. On the other hand, 


176 Differential Equations in Applications 


a vortex point is a fine singular point; its 
nature changes even under small perturba- 
tions of the right-hand sides. 


2.0 Stability of Equilibrium Points 
and of Periodic Motion 


As we already know, singular points of 
different types are characterized by different 
patterns of the phase trajectories in sufficient- 
ly small neighborhoods of these points. There 
is also another characteristic, the stability 
of an equilibrium point, which provides 
additional information on the behavior of 
phase trajectories in the neighborhood of 
singular points. Consider the pendulum de- 
picted in Figure 57. Two states of equilib- 
rium are shown: (a) an object of mass m is 
in a state of equilibrium at the uppermost 
point, and (b) the same object is in a state 
of equilibrium at its lowest point. The 
first state is unstable and the second, stable. 
And this is why. If the object is in its 
uppermost state of equilibrium, a slight 
push is enough to start it moving with an 
ever increasing speed away from the equilib- 
rium position and, hence, away from the 
initial position. But if the object is in the 
lowest possible state, a push makes it move 
away from the position of equilibrium 
with a decreasing speed, and the weaker the 
push the smaller the distance by which the 


= 


Ch. 2. Qualitative Methods 477 


Fig. 57 


object is displaced from the initial posi- 
tion. 

The state of equilibrium of a physical 
system corresponds to a singular point in 
the phase plane. Small perturbations at an 
unstable point of equilibrium lead to large 
displacements from this point, while at 
a stable point of equilibrium small pertur- 
bations lead to small displacements. Start- 
ing from these pictorial ideas, let us 
consider an isolated singular point of sys- 
tem (122), assuming for the sake of simplici- 
ty that the point is at the origin O (0, 0) 
of the phase plane. We will say that this 
singular point is stable if for every positive R 
there exists a positive r< R such that 
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Nikolai G. Chetaev (1902-1959) wrote: * 


...1f a passenger plane is being designed, a 
certain degree of stability must be provided 
for in the future movements of the plane so 
that it will be stable in flight and accident- 
free during take-off and landing. The crank- 
shaft must be so designed that it does not 
break from the vibrations that appear in real 
conditions of motor operation. To ensure that 
an artillery gun has the highest possible ac- 
curacy of aim and the smallest possible spread, 
the gun and the projectiles must be con- 
structed in such a manner that the trajectories 
of projectile flight are stable and the projec- 
tiles fly correctly. 

Numerous examples can be added to this 
list, and all will prove that real movements 
require selecting out of the possible solutions 
of the equations of motion only those that 
correspond to stable states. Moreover, if we 
wish to avoid a certain solution, it is advisable 
to change the design of the corresponding 
device in such a way that the state of motion 
corresponding to this solution becomes un- 
stable. 


Returning to the pendulum depicted in 
Figure 57, we note the following curious 
and somewhat unexpected fact. Research 
has shown that the upper (unstable) posi- 
tion of equilibrium can be made stable by 


* See N.G. Chetaev, Stability of Motion (Moscow: 
Nauka, 1965: pp. 8-9 (in Russian). 
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introducing vertical oscillations of the 
point of suspension. More than that, 
not only the upper (vertical) position of 
equilibrium can be made stable but also any 
other of the pendulum’s positions (for one, 
a horizontal position) by properly vibrating 
the point of suspension. * 

Now let us turn to a concept no less im- 
portant than the stability of an equilibrium 
point, the concept of stability of periodic 
movements (solutions). Let us assume that 
we are studying a conservative system that 
has periodic solutions. In the phase plane 
these solutions are represented by closed 
trajectories that completely fill a certain 
region. Thus, to each periodic motion of 
a conservative system there corresponds 
a motion of the representative point along 
a fixed closed trajectory in the phase plane. 

Generally, the period of traversal of 
different trajectories by representative 
points is different. In other words, the period 
of oscillations in a conservative system 
depends on the initial data. Geometrically 
this means that two closely spaced repre- 
sentative points that begin moving at 
a certain moment ¢ = f, (say, at the x axis) 
will move apart to a certain finite distance 


* The reader can find many cxarrplcs cf stakili- 
zation of differcnt ty pcs of pcndulums in the book 
by T.G. Strizhak, Methods of Investigating Dynam- 
ical Systems of the “Pendulum” Type (Alma-Ata: 
Nauka, 1981) (in Russian). 
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Fig. 58 


every phase trajectory originating at the 
initial moment ¢t = ¢, at a point P lying 
inside the circle 2? + y*? =r? will lie 
inside the circle zx? + y? = R?* for all t > fy 
(Figure 58). Without adhering to rigorous 
reasoning we can say that a singular point 
is stable if all phase trajectories that are 
near the point initially remain there with 
the passage of time. Also, a singular point 
is said to be asymptotically stable if it is 
stable and if there exists a circle 2? + 
y?> =r; such that each trajectory that at 
time t = ft, lies inside the circle converges 
to the origin as t+ -+oo, Finally, if 
a singular point is not stable, it is said to 
be unstable. 
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A singular point of the vortex type is 
always stable (but not asymptotically). 
A saddle point is always unstable. In 
Figure 59, which illustrates the behavior 
of the phase trajectories for pendulum oscil- 
lations in a medium with low drag, the sin- 
gular points, focal points, are asymptotical- 
ly stable; in Figure 56 the singular points, 
which are nodal points, are also asymptoti- 
cally stable. 

The introduced concept of stability of an 
equilibrium point is purely qualitative, 
since no mention of properties referring to 
the behavior of phase trajectories has been 
made. As for the concept of asymptotic 
stability, if compared with the notion of 
simple stability, it is additionally neces- 
sary that every phase trajectory tend to 
the origin with the passage of time. How- 
ever, in this case, too, no conditions are 
imposed on how the phase trajectory must 
approach point O (0, 0). 

The concepts of stability and asymptotic 
stability play an important role in appli- 
cations. The fact is that if a device is 
designed without due regard for stability 
considerations, when built it will be sensi- 
tive to the very smallest external perturba- 
tions, which in the final analysis may lead 
to extremely unpleasant consequences. Em- 
phasizing the importance of the concept of 
stability, the well-known Soviet specialist 
in the field of mathematics and mechanics 
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Fig. 59 


with the passage of time. But it also may 
happen that these points do not separate. 
To distinguish between these two possibil- 
ities, the concept of stability in the sense 
of Lyapunov is introduced for periodic 
solutions. The essence of this concept lies 
in the following. If knowing an e-neighbor- 
hood (with ¢ as smal] as desired) of a point 
M moving along a closed trajectory I 
(Figure 59) * ensures that we: know a mov- 
ing 6 (e)-neighborhood of the same point M@ 
such that every representative point that 
initially lies in the 6 (e)-neighborhood will 
never leave the e-neighborhood with the 
passage of time, then the periodic solution 


* An e-neighborhood of a point M is understood 
to be a disk of radius ¢ centered at point MM, 
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corresponding to [ is said to be stable in 
the sense of Lyapunov. If a periodic solution 
is not stable in the sense of Lyapunov, it is 
said to be unstable in the sense of Lyapunov. 

When it comes to periodic solutions that 
are unstable in the sense of Lyapunov, we 
must bear in mind that they still possess 
some sort of stability, orbital stability, 
which means that under small variations of 
initial data the representative point trans- 
fers from one phase trajectory to another 
lying as close as desired to the initially 
considered trajectory. 

Examples of periodic solutions that are 
stable in the sense of Lyapunov are those 
that emerge, for instance, when we consider 
the differential equation that describes the 
horizontal movements of an object of mass 
m in a vacuum with two linear springs acting 
on the object (Figure 52). An example of 
periodic solutions that are unstable in the 
sense of Lyapunov but are orbitally stable 
is the solutions of the differential equa- 
tion (134), which describes the motion of 
a circular pendulum in a medium without 
drag. 

In the first case the oscillation period 
does not depend on the initial data and is 
found by using the formula 7 = 2nV m/k. 
In the second the oscillation period depends 
on the initial data and is expressed, as we 
know, in terms of an elliptic integral of the 
first kind taken from Q to 3/2, 


484 Differential Equations in Applications 


Finally, we note that the question of 
whether periodic movements are stable in 
the sense of Lyapunov is directly linked 
to the question of isochronous vibrations. * 


2.6 Lyapunov Functions 


Intuitively it is clear that if the total 
energy of a physical system is at its mini- 
mum at a point of equilibrium, the point 
is one of stable equilibrium. This idea lies 
at the base of one of two methods used in 
studying stability problems, both suggested 
by the famous Russian mathematician Alek- 
sandr M. Lyapunov (1857-1918). This meth- 
od is known as Lyapunov’s direct, or 
second, method for stability investigations. ** 

We illustrate Lyapunov’s direct method 
using the (122) type of system when the 
origin is a singular point. 

Suppose that I’ is a phase trajectory of 
system (122). We consider a function V = 
V (zx, y) that is continuous together with its 
first partial derivatives 0V/dx and dV/dy 


* See, for example, the book by V. V. Amel’kin, 
N. A. Lukashevich, and A.P. Sadovskii, Non- 
linear Vibrations in Second-Order Systems (Minsk: 
Belorussian Univ. Press, 1982) (in Russian). 

** The reader will find many interesting examples 
of stability investigations involving differential 
models in the book by N. Rouche, P. Habets, and 
M. Laloy, Stability Theory by Lyapunov'’s Direct 
Method (New York: Springer, 1977). 
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in a domain containing I in the phase plane. 
If the representative point (z (ft), y (t)) 
moves along curve IT, then on this pee 
the function V (z, y) may be considered 
a function of ¢t, with the result that the 
rate at which V (zx, y) varies along I is 
given by the formula 


dV ov d oV d 
Gea te Get ty at ae XY) 
av 
Ta, Y(t y); (136) 


where X (zx, y) and Y (za, y) are the right- 
hand sides of system (422). 

Formula (136) is essential in the realiza- 
tion of Lyapunov’s direct method. The 
following concepts are important for the 
practical application of the method. 

Suppose that V = V (z, y) is continu- 
ous together with its first partial deriva- 
tives OV/dx and OV/dy in a domain G con- 
taining the origin in the phase plane, with 
V (0, 0) = 0. This function is said to be 
positive (negative) definite if at all points 
of G except the origin V (z, y) it is positive 
(negative). But if at points of G we have 
V (x, y) > 0 (< 0), the function V = V (2, 
y) is said to be nonnegative (nonpositive). 
For example, the function V defined by the 
formula V (xz, y) = x* + y® and considered 
in the (z, y)-plane is positive definite, while 
the function V (z, y) = x? is nonnegative 
since it vanishes on the entire y axis, 
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If V (x, y) is positive definite, we can 
require that at all points of the domain 
GNO the following inequality hold true: 


OV \2 OV \2 

(a2) +(4,) #2 

This means that the equation z = V (a, y) 
can be interpreted as the equation of a sur- 
face resembling a paraboloid that touches 
the (z, y)-plane at point O (0,0) (Figure 
60). Generally, the equation z = V (a, y 
with a positive definite V may specify 
a surface of a more complex structure. One 
such surface is shown in Figure 61, where 
the section of the surface with the plane 
z =C results not in a curve but in a ring. 


If a positive definite function V (2, y) is 
such that 


W (x, y)= ee 2 OX (2, y) 
OV (zx, 
ere PY (x, y) (137) 


is nonpositive, V is said to be the Lyapunov 
function of system (122). We note here that 
in view of (136) the requirement that W be 
nonpositive means that dV/dt< 0 and, 
hence, the function V = V (z, y) does not 
increase along the trajectory [ in the 
neighborhood of the origin. 

Here is a result arrived at by Lyapunov: 
iffor system (122) there exists a Lyapunov 
junction V (zx, y), then the origin, which is 
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Fig. 60 


Fig. 61 


a singular point, is stable. If the positive 
definite function V = V (az, y) is such that 
function W defined via (137) is negative defi- 
nite, then the origin is asymptotically stable. 

We will show with an example how to 
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apply the above result. Let us consider the 
equation of motion of an object of unit 
mass under a force exerted by a spring, 
which in view of (124) can be written in the 
form 
d2z 


as +c = 7 +kr=0, c>0. (138) 


The reader will recall that in this equation 
c > 0 characterizes the drag of the medium 
in which the object moves and k > 0 
characterizes the properties of the spring 
(the spring constant). The autonomous sys- 
tem corresponding to Eq. (138) has the 
form 


d 
Say, Ya —ke—cy. (139) 


For this system the origin of the phase 
(z, y)-plane is the only singular point. The 
kinetic energy of the object (of unit mass) 
is y?/2 and the potential energy (i.e. the 
energy stored by the spring) is 


( aeag = hee 


0 


This implies that the total energy of the 
system is 


V (x, Y= y ys kz. (140) 


It is easy to see that V specified by (140) 
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is a positive definite function. And since 
in the given case 


aV ov 
oe X(t YW) + ZY (a, y) 
=kzry +y(—kxr—cy) = —cy? <0, 


this function is the Lyapunov function of 
system (139), which means that the singu- 
lar point O (0, 0) is stable. 

In the above example the result was 
obtained fairly quickly. This is not always 
the case, however. The formulated Lyapu- 
nov criterion is purely qualitative and 
this does not provide a procedure for 
finding the Lyapunov function even if we 
know that such a function exists. This 
makes it much more difficult to determine 
whether a concrete system is stable or not. 

The reader must bear in mind that the 
above criterion of Lyapunov must be seen 
as a device for finding effective indications 
of equilibrium. Many studies have been 
devoted to this problem and a number of 
interesting results have been obtained in 
recent years.* 


2.7 Simple States of Equilibrium 


The dynamical interpretation of second- 
order differential equations already implies 
* The interested reader can refer to the book by 


E.A. Barbashin, Lyapunov Functions (Moscow: 
Nauka, 1970) (in Russian). 
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that investigation of the nature of equilibri- 
um states or, which is the same, the singular 
points provides a key for establishing the 
behavior of integral curves. 

It is also clear that, generally speaking, 
differential equations cannot be integrated 
in closed form. We need criteria that will 
enable us to determine the type of a sin- 
gular point from the form of the initial 
differential equation. Unfortunately, as 
a rule it is extremely difficult to find such 
criteria, but it is possible to isolate certain 
classes of differential equations for which 
this can be done fairly easily. Below we 
show, using the example of an object of unit 
mass subjected to the action of linear 
springs and moving in a medium with 
linear drag, how some results of the quali- 
tative theory of differential equations can 
be used to this end. But first let us discuss 
a system of the (122) type. It so happens 
that the simplest case in establishing the 
type of asingular point is when the Jacobian 
or functional determinant 


8X ox 
Ox Oy 
T(t, Y=! ay ay 
Oz Oy 


is nonzero at the point. 

If (z*, y*) isasingular point of system (4122) 
and if J* = J (z*, y*) «0, then the type 
of the singular point, which in the given 
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case is called a simple singular point, 
depends largely on the sign of constant J*. 
For instance, if J/* is negative, the singular 
point («*, y*) is a saddle point, and if J* is 
positive, the singular point may be a vor- 
tex point, a nodal point, or a focal point. 
The singular point may be a vortex point 
only if the divergence 


OX oY 
Dit, WHAT a, 


vanishes at the singular point, that is, only 
if D* = D (x*, y*) = 0. Note, however, 
that the condition D* =O is generally 
insufficient for the singular point (z*, y*) to 
be a vortex point. For a vortex point to 
be present certain additional conditions 
must be met, conditions that include high- 
er partial derivatives. And, generally, 
there can be an infinite number of such 
conditions. But if functions X and Y are 
linear in variables xz and y, the condition 
D* = 0 becomes sufficient for the singular 
point (x*, y*) to be a vortex point. 

If J/* > 0 but the singular point (z*, y*) 
is not a vortex point, a sufficiently small 
neighborhood of this point is filled with 
trajectories that either spiral into this point 
or converge to it in certain directions. Here, 
if D* > 0, the singular point is reached as 
t—» —oo and is unstable while if D* <0, 
the singular point is reached as t-—» -+ 00 
and proves to be stable. If the phase tra- 


192 Differential Equations in Applications 


jectories that reach the singular point are 
spirals, we are dealing with a focal point, 
but if the integral curves converge to a sin- 
gular point along a certain tangent, the 
point is a nodal point (Figure 62). 

Irrespective of the sign of Jacobian J*, 
the tangents to the trajectories of the differ- 
ential system (122) at a singular point (2*, 
y*) can be found from what is known as the 
characteristic equation 


a¥ (e*) yt) ~ | AY (et, yt) ~ 


y a Ox ~ Oy 

— ox (x*, y*) ed ax (z*, y*) ~ ? (141) 
Ox zp Oy y 

where 

r= 2 —2*, y=y—y*. (142) 


If X and Y contain linear terms, the 
partial derivatives in Eq. (141) act as co- 
efficients of z and y in the system obtained 
from system (122) after introducing the 
substitution (442). 

Equation (141) is homogeneous. Hence, 


if we introduce the slope A = yle of what is 
known as exceptional directions, we have 
the following quadratic equation for find- 
ing X: 

XM + (XZ —Ye)A— YE = 0. (143) 
The discriminant of this equation is 


A = (X2-+ Y¥%)? — 4J* = D*? — 4J*, 


13-0770 
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Hence, if J* <0, which corresponds to 
a saddle point, Eq. (143) always gives two 
real exceptional directions. But if J* > 0, 
there are no real exceptional directions, 
or there are two, or there is one 2-fold 
direction. In the first case the singular 
point is either a vortex point or a focal point. 

The existence of real exceptional direc- 
tions means (provided that J* is positive) 
that there is a singular point of the nodal 
type. For one thing, if there are two real 
exceptional directions, it can be proved 
that there are exactly two trajectories (one 
oneach side) whose tangent at the singular 
point is one of the exceptional straight lines 
(directions) while all the other trajecto- 
ries “enter” the singular point touching the 
other exceptional straight line (Figure 62a). 

If A = 0 and Eq. (143) is not an identity, 
we have only one exceptional straight line. 
The pattern of the trajectories for this case 
is illustrated by Figure 62b. It can be 
obtained from the previous case when the 
two exceptional directions coincide. The 
singular point divides the exceptional 
straight line into two half-lines, J, and l,, 
while the neighborhood of the singular point 
is divided into two sectors, one of which is 
completely filled with trajectories that 
“enter” the singular point and touch J, and 
the other is completely filled with trajecto- 
ries that “enter” the singular point and 
touch Z,. The boundary between the sectors 


= 
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consists of two trajectories, one of which 


‘toudnes’s, au tie Simputarporivana tne étner 


touches 1, at this point. 

If in Eq. (143) all coefficients vanish, we 
arrive at an identity, and then all the 
straight lines passing through the singular 
point are exceptional and there are exactly 
two trajectories (one on each side) that 
touch each of these straight lines at the 
singular point. This point (Figure 62c) is 
similar to the point with one 2-fold real 
exceptional direction. 


2.8 Motion of a Unit-Mass Object 
Under the Action of Linear Springs 
in a Medium with Linear Drag 


As demonstrated earlier, the differential 
equation describing the motion of a unit- 
mass object under the action of linear 
springs in a medium with linear drag has 
the form 


Se tes + ke=0. (144) 


So as not to restrict the differential model 
(144) to particular cases we will] not fix the 
directions in which the forces —c (dz/dt) 
and ~—kx act. As shown earlier, with 
Eq. (144) we can associate an autonomous 
system of the form 

dz d 

aM Gr heey. (145) 
13% 
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If we now exclude the trivial case with 
k =0, which we assume is true for the 
meantime, the differential system (145) has 
an isolated singular point at the origin. 
System (145) is a particular case of the 
general system (122). In our concrete exam- 
ple 


X (x, y) =y, Y (x, y) = —kr — cy, 
the Jacobian J (z, y) =k, and the divergence 


D (xz, y) = —c. The characteristic equa- 
tion assumes the form 
M+t+echtkh=0O, 


where A = c® — 4k is the discriminant of 
this equation. In accordance with the results 
obtained in Section 2.7 we arrive at the 
following cases. 

(1) If & is negative, the singular point is 
a saddle point with one positive and one 
negative exceptional direction. The phase 
trajectory pattern is illustrated in Fig- 
ure 63, where we can distinguish between 
three different types of motion. When the 
initial conditions correspond to point a, 
at which the velocity vector is directed to 
the origin and the velocity is sufficiently 
great, the representative point moves along 
a trajectory toward the singular point at 
a decreasing speed; after passing the origin 
the representative point moves away from 
it at an increasing speed. If the initial 
velocity decreases to a critical value, which 


= 
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corresponds to point b, the representative 
point approaches the singular point at 
a decreasing speed and “reaches” the origin 
in an infinitely long time interval. Finally, 
if the initial velocity is lower than the 
critical value and corresponds, say, to point 
C, the representative point approaches the 
origin at a decreasing speed, which vanishes 
at a certain distance z, from the origin. 
At point (z,, 0) the velocity vector reverses 
its direction and the representative point 
moves away from the origin. 

If the phase point corresponding to the 
initial state of the dynamical system lies 
in either one of the other three quadrants, 
the interpretation of the motion is obvious. 
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lig. 64 


Fig. 65 


(2) If k > 0, then J* is positive and the 
type of singular point depends on the 
value of c. This leaves us with the fol- 
lowing possibilities: 

(2a) If c = 0, that is, drag is nil, the 
singular point is a vortex point (Figure 64). 
The movements are periodic and their 
amplitude depends on the initial condi- 
tions, 
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(2b) Ifc > 0, that is, damping is posi- 
tive, the divergence D = 0X/dx + OY /dy is 
negative and, hence, the representative 
point moves along a trajectory toward the 
origin and reaches it in an infinitely long 
time interval. 

More precisely: 

(2b,) If A <0, that is, c? < 4k, the sin- 
gular point proves to be a focal point 
(Figure 65) and, hence, the dynamical 
system performs damped oscillations about 
the state of equilibrium with a decaying 
amplitude. 

(2b,) If A=O, that is, c? = 4k, the 
singular point is a nodal point witha single 
negative exceptional direction (Figure 66). 
The motion in this case is aperiodic and cor- 
responds to the so-called critical damping. 

(2b,) If A>O, that is, c> > 4k, the 
singular point is a nodal point with two 
negative exceptional directions (Figure 67). 
Qualitatively the motion of the dynamical 
system is the same as in the previous case 
and corresponds to damped oscillations. 

From the above results it follows that 
when c >0 and k>0, that is, drag is pos- 
itive and the restoring force is attractive, 
the dynamical system tends to a state of 
equilibrium and its motion is stable. 

(2c) Ife << 0, that is, damping is nega- 
tive, the qualitative pattern of the phase 
trajectories is the same as in the case (2b), 
the only difference being that here the 
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Fig. 66 


Fig. 67 


= 
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Fig. 69 


dynamical system ceases to be stable. 
Figure 68 contains all the above results 
and the dependence of the type of singular 
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point on the values of parameters c and k. 
Note that the diagrams can also be inter- 
preted as a summary of the results of studies 
of the types of singular points of system (122) 


When "9° =SSu-ac = “Hr aha’ « ="5r. 
However, the fact that c =0O does not 
generally mean that system (122) possesses 
a vortex point, and the fact that k = 0 does 
not mean that the system of a general type 
has no singular point. These cases belong 
to those of complex singular points, which 
we consider below. 

Returning to the dynamical system con- 
sidered in this section, we note that if k =0 
(c 0), the autonomous system (145) cor- 
responding to Eq. (144) assumes the form 


Sey Wes ay 
deo! ar ee 


This implies that the straight line y = 0 is 
densely populated by singular points, with 
the phase trajectory pattern shown in 
Figure 69. 

Finally, ifk =c = 0, Eq. (144) assumes 
the form 


The respective autonomous system is 


dz dy 
ae) ap 
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Here, as in the previous case, the z axis is 
densely populated by singular points. The 
respective phase trajectory pattern is 
shown in Figure 70. 


2.9 Adiabatic Flow 
of a Perfect Gas Through a Nozzle 
of Varying Cross Section 


A study of the flow of compressible viscous 
media is highly important from the practi- 
cal viewpoint. For one thing, such flow 
emerges in the vicinity of a wing and fuse- 
lage of an airplane; it also influences the 
operation of steam and gas turbines, jet 
engines, the nuclear reactors. 

Below we discuss the flow of a perfect 
gas through a nozzle with a varying cross 
section (Figure 71); the specific heat capaci- 
ty of the gas is cp and the nozzle’s varying 
cross-sectional area is denoted by A. The 
flow is interpreted as one-dimensional, that 
is, all its properties are assumed to be uni- 
form in a single cross section of the nozzle. 
Friction in the boundary Jayer is caused 
by the tangential stress t given by the for- 
mula 


Tt = gov"/2, (146) 
where g is the friction coefficient depending 


basically on the Reynolds number but 
assumed constant along the nozzle, p the 
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Y=t 


Fig. 70 
=: 


Fig. 71 


flux density, and v the flow velocity. Final- 
ly, we assume that adiabaticity conditions 
are met, that is, drag, combustion, chemi- 
cal transformations, evaporation, and con- 
densation are excluded. 

One of the basic equations describing 
this type of flow is the well-known conti- 
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nuity equation, which in the given case 
is written in the form 


w = pAv, (147) 


where the flux variation rate wis assumed 
constant. From this equation it follows 
that 


dp dA dv 
et ae a (148) 

Let us now turn to the equation for the 
energy of steady-state flow. We note that 
generally such an equation links the exter- 
nal work done on the system and the action 
of external heat sources with the increase 
in enthalpy (heat content) in the flow and 
the kinetic and potential energies. In our 
case the flow is adiabatic; hence, the energy- 
balance equation can be written in the form 


O = w (h + dh) — wh 

+ w [v?/2 +d (v*/2)] — wv*/2, 

or 

dh + d (v?/2) = 0, (149) 


where h is the enthalpy of the flow (the 
thermodynamic potential) at absolute tem- 
perature 7. But in Eq. (149) dh =c,dT 
and, therefore, we can write the equation 
for the flow energy as 


c,dT + d (v/2) = 0. (150) 
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Now let us derive the momentum equa- 
tion for the flow. Note that here the common 
approach to problems involving a steady- 
state flow is to use Newton’s second law. 
Assuming that the divergence angle of the 
nozzle wall is small, we can write the 
momentum equation in the form 


pA + pdA — (p + dp) (A + dA) — tdA 


= w dv, 
or 
—A dp—dAdp—tdA =wd, (151) 


where p is the static pressure. 

The term dA dp in Eq. (151) is of a 
higher order than the other terms and, 
therefore, we can always assume that the 
momentum equation for the flow has the 
form 


—A dp —tdA =—w dv. (152) 


If we denote by D the hydraulic diameter, 
we note that its variation along the nozzle’s 
axis is determined by a function F such that 
D = F (x), where z is the coordinate along 
the nozzle’s axis. From the definition of the 
hydraulic diameter it follows that 


dA 4dz 
A D° 


(153) 


Noting that pv?/2 = ypM?*, where y is the 
specific heat ratio of the medium, and M 
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is the Mach number, we can write formula 
(146) thus: 


t = gypM?. (154) 


Combining this with Eqs. (147) and (153), 
we arrive at the following represen- 
tation for the momentum equation (152): 


dp dz dv2 


Denoting the square of the Mach number by 
y, employing Eqs. (148), (150), and (155), 
and performing the necessary algebraic 
transformations, we arrive at the following 
differential equation: 


ay fu (442) oF) a 
dz (i—y) F (a) yee) 


where the prime stands for derivative of 
the respective quantity.* 

The denominator of the right-hand side of 
Eq. (156) vanishes at y = 1, that is, when 
the Mach number becomes equal to unity. 
This means that the integral curves of the 
last equation intersect what is known as the 
sonic line and have vertical tangents at the 
intersection points. Since the right-hand 


* See J. Kestin and S.K. Zaremba, “One-dimen- 
sional high-speed flows. Flow patterns derived 
for the flow of gases through nozzles, including 
compressibility and viscosity effects,” Aircraft 
Engin. 25, No. 292: 172-175, 179 (1953). 
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side of Eq. (156) reverses its sign in the 
process of intersection, the integral curves 
“flip over” and the possibility of inflection 
points is excluded. The physical meaning 
of this phenomenon implies that along 
integral curves the value of x must increase 
continuously. Hence, the section on which 
the integral curves intersect the sonic line 
with vertical tangents must be the exit 
section of the nozzle. Thus, the transition 
from subsonic flow to supersonic (and back) 
can occur inside the nozzle only through 
a singular point with real exceptional di- 
rections, that is, through a saddle point or 
a nodal point. 

The coordinates of the singular points of 
Eq. (196) are specified by the equations 


y* =1, F (x*) = 9, 


which imply that these points are situated 
in the diverging part of the nozzle. A saddle 
point appears if /* is negative, that is, 
F" (z*) > 0. Since gq is a sufficiently small 
constant, a saddle point appears near the 
throat of the nozzle. A nodal point, on the 
other hand, appears only if F”"(z*) is neg- 
ative. Thus, a nodal point emerges in the 
part of the nozzle that lies behind an in- 
flection point of the nozzle’s profile or, in 
practical terms, at a certain distance from 
the throat of the nozzle, provided that the 
profile contains an inflection point. 


Ch. 2. Qualitative Methods 209 


From the characteristic equation 
F (x*) h* + 2gqy (y + 1) A 
—2(y +1) F"(c*) =0 


we can see that the slopes of the two excep- 
tional directions are opposite in sign in the 
case of a saddle point and have the same 
sign (are negative) in the case of a nodal 
point. This means that only a saddle point 
allows for a transition from supersonic to 
subsonic velocities and from subsonic to 
supersonic velocities (Figure 72). The case 
of a nodal point (Figure 73) allows for a 
continuous transition only from supersonic 
to subsonic flow. 

Since Eq. (156) cannot be integrated in 
closed form, we must employ numerical 
methods of integration in any further dis- 
cussion. It is advisable in this connection 
to begin the construction of the four sepa- 
ratrices of a saddle point as integral curves 
by allowing for the fact that the singular 
point itself is a point, so to say, from which 
these integral curves emerge. Such a con- 
struction is indeed possible since the char- 
acteristic equation provides us with the 
direction of the two tangents at the singular 
point S (x*, 1). If this is ignored and the 
motion monitored as beginning at points 
a and b in Figure 72, which lie on different 
sides of a separatrix, the corresponding 
points move along curves a and f that 
strongly diverge and, hence, provide no in- 
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Fig. 73 


formation about the integral curve (the 
separatrix) that “enters” point S. On the 
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other hand, if we move along an integral 
curve that “emerges” from point S and we 
assume that the initial segment of the curve 
coincides with a segment of an exceptional 
straight line, the error may be minimized 
if we allow for the convergence of integral 
curves in the direction in which the values 
of x diminish. 

Figure 72 illustrates the pattern of inte- 
gral curves in the vicinity of a singular 
point. The straight line passing through 
point xz; (the throat of the nozzle) corre- 
sponds to values at which the numerator 
in the right-hand side of Eq. (156) vanishes, 
which points to the presence of extrema. 


2.10 Higher-Order Points 
of Equilibrium 


In previous sections we studied the types 
of singular points that emerge when the 
Jacobian J* is nonzero. But suppose that 
all the partial derivatives of the functions X 
and Y in the right-hand sides of system (122) 
vanish up to the mth order inclusive. Then 
in the vicinity of a singular point there 
may be an infinitude of phase-trajectory 
patterns. However, if we exclude the pos- 
sibility of equilibrium points of the vortex 
and focal types emerging in this picture, 
then it appears that the neighborhood of 
a singular point can be broken down into 


14* 
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a finite number of sectors belonging to 
three standard types. These are the hyper- 
bolic, parabolic, and elliptic. Below we de- 
scribe these sectors in detail, but first let us 
make several assumptions to simplify mat- 
ters. 

We assume that the origin is shifted to 
the singular point, that is, z* = y* = 0; 
the right-hand sides of system (122) can 
be written in the form 


x (x, y) =X, (x, y) + @ (x, y), 
Y@Q.y=Yi@,y+¥@y, (9 


where X, and Y, are polynomials of de- 
gree n homogeneous in variables x and y 
(one of these polynomials may be identically 
zero), and the functions M and W have in 
the neighborhood of the origin continuous 
first partial derivatives. In addition, we 
assume that the functions 


@D (z, y) @, (z, y) D, (z, y) 
(22 y2)(M*H/2 % (at y2)"/2 7 “Ga anya 
¥ (z, y) W. (2, y) Vy (x, y) 


(x2 y2)(™*I)/2 7 (a2 2)? aya 


are bounded in the neighborhood of the 
origin. Under these assumptions the follow- 
ing assertions hold true. 

(1) Every trajectory of the system of 
equations (122) with right-hand sides of the 
(157) form that “enters” the origin along 
a certain tangent touches one of the ex- 
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ceptional straight lines specified by the 
equation 


rY , (t, y) — yXn (z, y) = 0. (158) 


Since the functions X, and Y, are homo- 
geneous, we can rewrite Eq. (158) as an 
equation for the slope A = y/x. Then the 
exceptional straight line is said to be sin- 
gular if 


X, (t, y) = Y, (@, y) = 0 


on this straight line. Some examples of 
such straight lines are shown in Figure 62. 
The straight lines defined by Eq. (158) but 
not singular are said to be regular. 

(2) The pattern of the phase trajectories 
of (422) in the vicinity of one of the two 
rays that “emerge” from the origin and 
together form an exceptional straight line 
can be studied by considering a small disk 
(centered at the origin) from which we 
select a sector limited by two radii lying 
sufficiently close to the ray on both sides 
of it. Such a sector is commonly known as 
a standard domain. 

More than that, in the case of a regular 
exceptional straight line, which  corre- 
sponds to a linear factor in Eq. (108), the 
standard domain considered in a disk of 
a sufficiently small radius belongs to either 
one of two types: attractive or repulsive. 

(2a) The attractive standard domain is 
characterized by the fact that each tra- 
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Fig. 74 


Fig. 75 


jectory passing through it reaches the 
origin along the tangent that coincides 
with the exceptional straight line (Fig- 
ure 74). 

(2b) The repulsive standard domain is 
characterized by the fact that only one 
phase trajectory passing through it reaches 
the singular point along the tangent that 
coincides with the exceptional straight line. 
All other phase trajectories of (122) that 
enter the standard domain through the 
boundary of the disk leave the disk by 
crossing one of the radii that limit the 
domain (Figure 75). 
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Let us turn our attention to the following 
fact. If the disk centered at the origin is 
small enough, the two types of standard 
domains can be classified according to the 
behavior of the vector (X, Y) on the 
boundary of the domain. The behavior of 
vector (X, Y) can be identified here with 
the behavior of vector (X,, Y,). More 
than that, it can be demonstrated that if, 
aS assumed, a fixed exceptional direction 
does not correspond to a multiple root of 
the characteristic equation, the vector con- 
sidered on one of the radii that limit the 
domain is directed either inward or out- 
ward. Then if in the first case the vector 
considered on the part of the boundary 
of the standard domain that is the arc of 
the circle is also directed inward, and in 
the second case outward, the standard 
domain is attractive. But if the opposite 
situation is true, the standard domain is 
repulsive. It must be noted that in any case 
the vector considered on the part of the 
boundary that is the arc of the circle is 
always directed either inward or outward 
since it is almost parallel to the radius. 

Standard domains corresponding to sin- 
gular exceptional directions or multiple 
roots of the characteristic equation have 
a more complicated nature, but since to 
some extent they constitute a highly rare 
phenomenon, we will not describe them here. 

If we now turn, for example, to a saddle 
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point, we note that it allows for four 
repulsive standard domains. In the neigh- 
borhood of a singular point of the nodal 
type there are two attractive standard 
domains and two repulsive. 

(3) If there exist real-valued exceptional- 
direction straight lines, the neighborhood 
of a singular point can be divided into a 
finite number of sectors each of which is 
bounded by the two phase trajectories of 
(122) that “enter” the origin along definite 
tangents. Each of such sectors belongs to 
one of the following three types. 

(3a) The elliptic sector (Figure 76) con- 
tains an infinitude of phase trajectories 
in the form of loops passing through the 
origin and touching on each side of the 
boundary of the sector. 

(3b) The parabolic sector (Figure 77) 
is filled with phase trajectories that connect 
the singular point with the boundary of the 
neighborhood. 

(3c) The hyperbolic sector (Figure 78) 
is filled with phase trajectories that ap- 
proach the boundary of the neighborhood 
in both directions. 

More precisely: 

(4a) Elliptic sectors are formed between 
two phase trajectories belonging to two 
successive standard domains, both of which 
are attractive. 

(4b) Parabolic sectors are formed between 
two phase trajectories belonging to” two 
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Fig. 76 Fig. 77 


Fig. 78 


successive standard domains one of which 
is attractive and the other repulsive. All 
phase trajectories that pass through the 
Jatter domain touch at the singular point 
of the exceptional! straight line that defines 
the attractive domain. 
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(4c) Hyperbolic sectors are formed be- 
tween two phase trajectories belonging to 
two successive repulsive standard domains. 

For example, it is easy to distinguish the 
four hyperbolic sectors at a saddle point 
and the four parabolic sectors at a nodal 
point. Elliptic sectors do not appear in 
the case of simple singular points, where 
the Jacobian J* is nonzero. 

If a singular point does not allow for the 
existence of real exceptional directions, 
the phase trajectories in its neighborhood 
always possess the vortex or focal struc- 
ture.* 


2.141 Inversion with Respect 
to a Circle and Homogeneous 
Coordinates 


Above we described methods for establish- 
ing the local behavior of phase trajectories 
of differential systems of the (422) type in 
the neighborhood of singular points. And 
although in many cases all required infor- 
mation can be extracted by following these 
methods, there may be a need to study the 


* Methods that make it possible to distinguish 
between a vortex point and a focal point are dis- 
cussed, for example, in the book by V.V Amel’- 
kin, N.A. Lukashevich, and A.P. Sadovskii, 
Nonlinear Vibrations in Second-Order Systems 
(Minsk: Belorussian Univ. Press, 1982) (in Rus- 
sian). 
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trajectory behavior in infinitely distant 
parts of the phase plane, as x? + y? +> oo. 
A simple way of studying the asymptotic 
behavior of the phase trajectories of the 
differential system (122) is to introduce 
the point at infinity by transforming the 
initial differential system in an appropriate 
manner, say by inversion, which is defined 
by the following formulas: 


eran eee 
on ae By? 
ee ate 
(E= app’? I ape ). p22) 


Geometrically this transformation consti- 
tutes what has become known as inversion 
with respect to a circle and maps the origin 
into the point at infinity and vice versa. 
Transformation (159) maps every finite 
point M (zx, y) of the phase plane into point 
M’ (&, n) of the same plane, with points M 
and M’ lying on a single ray that emerges 
from the origin and obeying the condition 
OM x OM’ =? (Figure 79). It is well 
known that such a transformation maps 
circles into circles (straight lines are con- 
sidered circles passing through the point 
at infinity). For one thing, straight lines 
passing through the origin are invariant 
under transformation (159). Hence, the 
slopes of asymptotic directions are the 
slopes of tangents at the new origin € = 
yn =0. Note that in the majority of 
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Fig. 80 


cases the new origin serves as a singular 
point. The reasons for this are discussed 
below. 
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As for the question of how to construct 
in the old (x, y)-plane a curve that has a 
definite asymptotic direction, that is, has 
a definite tangent at the origin of the new 
(E, y)-plane, this can be started by con- 
sidering the (&, y)-plane, more precisely, 
by considering this curve in, say, a unit 
circle in the (&, y)-plane. The fact is that 
since a unit circle is mapped, via transfor- 
mation (159), into itself, we can always estab- 
lish the point where the curve intersects the 
respective unit circle in the (z, y)-plane; 
any further investigations can be carried 
out in the usual manner. 

We also note that completion of the 
(zr, y)-plane with the point at infinity is 
topologically equivalent to inversion of 
the stereographic projection (Figure 80), 
in which the points on a sphere are mapped 
onto a plane that is tangent to the sphere 
at point S. The projection center N is the 
antipodal point of S. It is clear that the 
projection center N corresponds to the 
point at infinity in the (x, y)-plane. Con- 
versely, if we map the plane onto the 
sphere, a vector field on the plane trans- 
forms into a vector field on the sphere and 
the point at infinity may prove to be a sin- 
gular point on the sphere. 

Although inversion with respect to a 
circle is useful, it proves cumbersome and 
inconvenient when the point at infinity has 
a complicated structure. In such cases 
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another, more convenient, transformation of 
the (z, y)-plane is used by introducing 
homogeneous coordinates: 


£2 Els. y=. 


Under this transformation, each point of 
the (z, y)-plane is associated with a triple 
of real numbers (€, yn, z) that are not 
simultaneously zero and no difference is 
made between the triples (&, y, z) and 
(KE, kn, kz) for every real k= 0. If a 
point (x, y) is not at infinity, z=4 0. But 
if z = 0, we have a straight line at infinity. 
The (x, y)-plane completed with the straight 
line at infinity is called the projective plane. 
Such a straight line may carry several sin- 
gular points, and the nature of these points 
is usually simpler than that of a singular 
point introduced by inversion with respect 
to a circle. 

If we now consider a pencil of lines and 
describe its center with a sphere of, say, 
unit radius, then each line of the pencil 
intersects the sphere at antipodal points. 
This implies that every point of the pro- 
jective plane is mapped continuously and 
in a one-to-one manner onto a pair of anti- 
podal points on the unit sphere. Thus, the 
projective plane may be interpreted as the 
set of all pairs of antipodal points on a unit 
sphere. To visualize the projective plane, 
we need only consider one-half of the sphere 
and assume its points to be the points of the 
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Fig. 84 


projective plane. If we orthogonally project 
this hemisphere (say, the lower one) onto 
the a-plane, which touches it at pole S 
(Figure 81), the projective plane is mapped 
onto a unit disk whose antipodal points 
on the boundary are assumed identical. 
Each pair of the antipodal points of the 
boundary corresponds to a line at infinity, 
and the completion of the Euclidean plane 
with this line transforms the plane into 
a closed surface, the projective plane. 


2.12 Flow of a Perfect Gas 
Through a Rotating Tube 
of Uniform Cross Section 


In some types of turboprop helicopters and 
airplanes and in jet turbines, the gaseous 
fuel-air mixture is forced through rotating 
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Fig. 82 


tubes of uniform cross section installed 
in the compressor blades and linked by a 
hollow vertical axle. To establish the opti- 
mal conditions for rotation we must analyze 
the flow of the mixture through a rotating 
tube and link the solution with boundary 
conditions determined by the tube design. 
In a blade the gaseous mixture participates 
in a rotational motion with respect to the 
axis with constant angular velocity 
and moves relative to the tube with an 
acceleration v (du/dr), where v is the speed 
of a gas particle with respect to the tube, 
and r is the coordinate measured along the 
rotating compressor blades. 

Figure 82 depicts schematically a single 
rotating tube of a compressor blade.* It is 


* See J. Kestin and S.K. Zaremba, “Adiabatic 
one-dimensional flow of a perfect gas through a 
rotating tube of uniform cross section,” Aeronaut. 
Quart. 4: 373-399 (1954). 
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assumed that the fuel-air mixture, whose 
initial state is known, is supplied along 
the hollow axle to a cavity on the axle in 
which the flow velocity may be assumed 
insignificant. The boundary of the cavity 
occupied by the gas is denoted by a, the 
gas is assumed perfect with a specific heat 
ratio y, and all the processes that the gas 
mixture undergoes are assumed reversibly 
adiabatic (exceptions are noted below). 

It is assumed that the gas expands through 
a nozzle with an outlet cross section b that 
at the same time is the inlet of a tube of 
uniform cross section. The expansion of 
the gas mixture from state a to state b 
is assumed to proceed isentropically; we 
denote the velocity of the gas after expan- 
sion by v, and the distance from the rotation 
axis O, to the cross section b by ry. 

When passing through the tube, whose 
uniform cross-sectional area will be denoted 
by A and whose hydraulic diameter by D, 
the gas is accelerated thanks to the com- 
bined action of the pressure drop and dy- 
namical acceleration in the rotation compres- 
sor blade. We ignore here the effect produced 
by pressure variations in the tube (if they 
exist at all) and by variations in pressure 
drop acting on the cross section plane, both 
of which are the result of the Coriolis force. 
The last assumption, generally speaking, 
requires experimental verification, since 
the existence of a lateral pressure drop may 
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serve as a cause of secondary flows. But if 
the diameter of the tube is small compared 
to the tube’s length, such an assumption 
is justified. 

It is now clear that the equations of mo- 
mentum and energy balance for a compres- 
sible mixture traveling along a tube of 
uniform cross section must be modified so 
as to allow for forces of inertia that appear 
in a rotating reference frame. As for the 
continuity equation, it remains the same. 

Now let us suppose that starting from 
cross section c at the right end of the tube 
the gas is compressed isentropically and 
passes through a diverging nozzle. In the 
process it transfers into a state of rest with 
respect to the compressor blade in the second 
cavity at a distance r, from the rotation 
axis and reaches a state with pressure P, 
and temperature 7. 

From the second cavity the gaseous mix- 
ture expands isentropically into a converg- 
ing or converging-diverging nozzle in such 
a manner that it leaves the cavity at right 
angles to the tube’s axis. This produces 
a thrust force caused by the presence of 
a torque. 

Below for the sake of simplicity we as- 
sume that the exit nozzle is a converging 
one and has an exit (throat) of area A*. 
Denoting the external (atmospheric) pres- 
sure by P,, we consider two modes of pas- 
sage of the mixture through the nozzle. 
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The first applies. to a situation in which 
the P,-to-P, ratio exceeds the critical 
value, or 


Po/Pq > (2/(y + 1))Vr-d. 


In this case the flow at the nozzle’s throat 
is subsonic and, hence, pressure P, at the 
throat is equal to the atmospheric pressure, 
that is, P, = P,. The second case applies 
to a situation in which the P,-to-P, ratio 
is below the critical value, or 


Pal Pa < (2/(y + 1)Vr-), 


The pressure P, at the nozzle’s throat has 
a fixed value that depends on P, but not 
on P,. Hence, 


P, = (2/(y + 1))/V-1) Pa. 


In the latter case the flow at the nozzle’s 
throat has the speed of sound 


vs = (2My + 1) aa, 


where a, depends only on temperature 7 j. 

In further analysis the flow is assumed 
adiabatic everywhere and isentropic every- 
where except in the tube between cross sec- 
tions 6 and c. 

As in the case where we derive the differ- 
ential equation that describes the adia- 
batic flow of a perfect gas through a nozzle 
of varying cross section, let us now consider 
the continuity equation, the equation of 
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momentum, and the equation of energy 
balance. For instance, the equation of con- 
tinuity in this case has the form 


pe ae const, (160) 


with V the specific volume and » the flux 
density, or 


» = m'/A, 


where m’ is the flux mass. 

To derive the equation of motion, or the 
momentum equation, we turn to Figure 83. 
We note that the dynamical effect of the 
rotational motion of the compressor blade 
can be used to describe the flow with re- 
spect to the moving tube, where in accord- 
ance with the D’Alembert principle the 
force of inertia 


_A 2 
dl = —- w*r dr 
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is assumed to act in the positive direction 
of r. Hence, the element of mass dm = 
(A/V) dr moves with an _ acceleration 
v (dv/dr) caused by the combined action 
of the force of inertia d/J, the force of pres- 
sure A dP, and the friction force dF = 
(4A/D) dr. Here t =A (v*?/2V), where A 
depends on the Reynolds number R. As 
a first approximation we can assume that 
X remains constant along the entire tube. 
With this in mind, we can write the mo- 
mentum equation in the form 


A dv ! =a 

7 drv-— = — AdP —=— dr r+ A wr dr, 
or after certain rN 

V dP +v dv+— vtdr—wrdr=0. (164) 


As for the energy-balance equation, it is 
simple to derive if we use the first law of 
thermodynamics for open systems and bear 
in mind that the amount of work performed 
by the system is w*r dr. Thus 


dh + v dv — w*r dr = 0, 


where h is the enthalpy. If we define the 
speed of sound a via the equation 


h= 


a2 


i oe 
we find that 
a da-+-— *vdv— 2—* wer dr =0. 
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This leads us to the following energy-balance 

equation: 

a2 = a2 PoE pe VE ape, (162) 

If we now introduce the dimensionless 

quantities 

M,=vlay, x =r/D, G = w*D?*/a?, (163) 

we arrive at the equation 

ae y—1 y—1 

Fr et mie eet Ce = 2. (164) 
By excluding the pressure and specific 

volume from the basic equations (4160), 

(161), and (164) we can derive an equation 

that links the dimensionless quantity M, 

with the dimensionless distance. This equa- 

tion serves as the basic equation for solving 

our problem. The continuity equation (160) 

implies that 


V = o/b. (165) 
Then, combining 


$i PY Gs y~Pv 
a y ~ 


with Eq. (162), we arrive at the following 
relationships: 


2 aa) — 2 
P= (Pvt ) yp, 


“yp 2y 2y ov 
dP y—1 wr {v= ! dis ax 
dr —sy v \ 29 yu 


p2 dr ° 


—4 wr2\ d a 
lg ae) (166) 
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Substituting the value of V from Eq. (165) 
and the value of dP/dr from Eq. (166) into 
the momentum equation (164) and allowing 
for the relationship (163), we arrive at the 
differential equation 


dMz 2M32 (2\yM2— G22) 


. (167 

dx ¢—VE1 yey Y—1 @az2 00 
2 2 
Assuming that 
1 
My=y, m=2hy, p=+(y+1), 
1 

we can rewrite Eq. (167) as 
dy _—«- 2y (my — Gz) 
de 1 py ta aed 


The differential equation (168) is the 
object of our further investigation. In- 
troducing the variables (159), we reduce it to 


dy 2m (22? — pn -+ gGE?)-+ 2y (n? — £2) Q 
dé ss A (S889 — py -} gGE?) ++ 46n?02 
(469) 


where 2 = €?-+ y7, Q= myn— GE, A= E? — x. 
Clearly, the origin constitutes a singular 
point for Eq. (169). The lowest-order terms 
in both numerator and denominator are 
quadratic. If we discard higher-order terms, 
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as we did on p. 211 when discussing higher- 
order equilibrium points, we find that 


Y, (8, n) 
= 2qG°Eq + 2n (yn? — §?) (mn — G8), 
X,4 (§, n) 
= qgG’E? (—? — n*) + 4En? (myn — GE), 


and, hence, the origin is a higher-order 
singular point. 

The characteristic equation of the differ- 
ential equation (169) assumes the following 
form after certain manipulations: 


En (6? + *) [(g + 2) GE — 2my] = 0. 
(170) 


This leads us to three real-valued excep- 
tional straight lines: 


— 0, b — 0, 
©) 2 ++ 2) ge 2my = 0. (171) 


Each is a regular straight line and corre- 
sponds to one of the factors in Eq. (170). 

It is easy to see that for & and y positive 
the value of X, (&, n) is positive in the 
neighborhood of all three exceptional 
straight lines, with the result that the 
expression 


m6 0 10 
X4 (, n) E 
__ En (E2-+ 4?) [(¢-+2) G2E—2my] (172) 


~———-& [G22 (E3 — 2) -+ 48n? (my — GE)] 
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has the same sign in both the first quadrant 
and the neighborhood of the straight lines 
as the left-hand side of (174c) and, hence, 
is negative below the straight line (174c) 
and positive above. The explanation lies 
in the fact that (172) may change sign only 
when an exceptional straight line is crossed. 
Geometrically this means that the right- 
hand side of the differential equation 


dy — Y, (&, n) 
d X4 (&, n) 


determining the behavior of integral curves 
in a disk of sufficiently small radius and 
centered at the origin fixes an angle greater 
than the angle of inclination of the radii 
that lie in the first quadrant between the 
straight lines (171b) and (1714c). As for 
the radii lying between the straight lines 
({71c) and (174a), on them the right-hand 
side of this differential equation fixes an 
angle that is smaller than the angle of 
inclination of these radii (Figure 84). Thus, 
the standard domain containing the straight 
line (474c) must be repulsive, while the 
domains containing the straight lines (171a) 
and ({471b) must be attractive. In view of 
the symmetry of the field specified by vector 
(X,, Y,), the above facts remain valid for 
standard domains obtained from those men- 
tioned earlier by a rotation about the 
singular point through an angle of 180°. 


234 Differential Equations in Applications 
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Thus, there are exactly two integral 
curves that “enter” the singular point along 
the tangent (171c) and an infinitude of 
integral curves that touch the coordinate 
axes ({71a) and (174b) at the point of rest. 

We see then that the second and fourth 
quadrants contain elliptic sectors since 
they lie between two successive attractive 
standard domains (Figure 85). Each of the 
first and third quadrants is divided into 
two sectors by the integral curves that 
touch the exceptional straight line (171c) 
at the origin. These sectors are parabolic 
because they lie between two successive 
standard domains, one of which is attractive 
and the other repulsive. More than that, 
all the integral curves except those that 
touch the straight line (171c) touch the 


= 
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coordinate axes (171a) and (471b) at the 
singular point. 

Basing our reasoning on physical con- 
siderations, we can analyze the differential 
equation (168) solely in the first quadrant 
of the (z, y)-plane. Turning to this plane, 
we see that there exists exactly one integral 
curve having an asymptotic direction with 
a slope (q + 2) G?/2m. All other integral 
curves allow for the asymptotic direction 
of one of the coordinate axes. Indeed, it is 
easy to prove that all these integral curves 
asymptotically approach one of the coor- 
dinate axes, that is, along each of them in 
the movement toward infinity not only does 
y/x +O or z/y +0 but so does y ~O or 
z—->Q (with +t-—->oo or y—-> oo, respec- 
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tively). This is illustrated in Figure 86, 
where some integral curves have been plot- 
ted at great values of z and y. Note that the 
pattern of integral curves at a finite dis- 
tance from the origin depends primarily 
on the position of the singular points, and 
a special investigation is required to clarify 
it. It can also be proved that the integral 
curves that approach the coordinate axes 
asymptotically (Figure 86) leave the first 
quadrant at a finite distance from the 
origin. 


2.13 Isolated Closed Trajectories 


We already know that in the case of a 
singular point of the vortex type a certain 
region of the phase plane is completely 
filled with closed trajectories. However, 
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a more complicated situation may occur 
when there is an isolated closed trajectory, 
that is, a trajectory in a certain neigh- 
borhood of which there are no other closed 
trajectories. This case is directly linked 
with the existence of isolated periodic 
solutions. Interestingly, only nonlinear dif- 
ferential equations and systems can have 
isolated closed trajectories. 

Isolated periodic solutions correspond 
to a broad spectrum of phenomena and pro- 
cesses occurring in biology, radiophysics, 
oscillation theory, astronomy, medicine, 
and the theory of device design. Such solu- 
tions emerge in differential models in eco- 
nomics, in various aspects of automatic 
control, in airplane design, and in other 
fields. Below we study the possibility of 
isolated periodic solutions emerging in pro- 
cesses that occur in electric circuits; we also 
consider as a model the nonlinear differ- 
ential system 


(173) 
dt ty —x—y’). 


To solve this system, we introduce polar 
coordinates r and 0, where x =rcos@ 
and y = rsin 0. Then, differentiating the 
relationships 2? +y?=r* and O= 
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tan-! (y/x) with respect to t, we arrive at 


dz d dr di dz dé 

SMe pls es Be OY ap BO 

Ce ee ap de te de ae 
(174) 


Multiplying the first equation into z and 
the second into y, adding the products, 
and allowing for the first relationship in 
(174), we find that 


r ST 72 (4 — 32). (175) 


Multiplying the second equation in (173) 
into z and the first into y, subtracting one 
product from another, and allowing for the 
second relationship in (174), we find that 


p? = 7?, (176) 


System (173) has only one singular point, 
O (0, 0). Since at the moment we are only 
interested in constructing trajectories, we 
can assume that r is positive. Then Eqs. 
(175) and (176) imply that system (173) 
can be reduced to the form 


A= 1. (177) 


Each of these equations can easily be inte- 
grated and the entire family of solutions, 
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as can easily be seen, is given by the for- 
mulas 


1 
VTL Gett ” 6=t+41,, (178) 
or, in terms of the old variables xz and y, 
_ _£0s (¢-+ to) = sin (t-+ fo) | 
V 1+ Ce-*t ’ V 14 Ce-2t 


If now in the first equation in (178) we put 
C=0, we get r=1 and 0=t+48. 
These two relationships define a closed 
trajectory, a circle z?+ y27 =1. If C is 
negative, it is clear that r is greater than 
unity and tends to unity as t ~ +00. But 
if C is positive, it is clear that r is less than 
unity and tends to unity as t — -+oo. This 
means that there exists only one closed 
trajectory r = 1 which all other trajectories 
approach along spirals with the passage of 
time (Figure 87). 

Closed phase trajectories possessing such 
properties are known as limit cycles or, 
more precisely, (orbitally) stable limit cycles. 
In fact, there can be two additional types 
of limit cycles. A limit cycle is said to be 
(orbitally) unstable if all neighboring tra- 
jectories spiral away from it as t — -+ oo. 
And a limit cycle is said to be (orbitally) 
half-stable if all neighboring trajectories 
on one side (say, from the inner side) spiral 
on to it and all neighboring trajectories 
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Fig. 87 


on the other side (say, on the outer side) 
spiral away from it as t +-boo. 

In the above example we were able to 
find in explicit form the equation of a closed 
phase trajectory, but generally, of course, 
this cannot be done. Hence the importanca 
in the theory of ordinary differential equa- 
tions of criteria that enable at least spe- 
cifying the regions where a limit cycle 
may occur. Note that a closed trajectory of 
(122), if such a trajectory exists, contains 
within its interior at least one singular 
point of the system. This, for one thing, 
implies that if there are no singular points 
of a differential system within a region of 
the phase plane, there are no closed tra- 
jectories in the region either. 


= 
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Let D be a bounded domain that lies 
together with its boundary in the phase 
plane and does not contain any singular 
points of system (122). Then the Poincaré- 
Bendixson criterion holds true, namely, if 
. is a trajectory of (122) that at the initial 
moment t = t, emerges from a point that 
lies in D and remains in D for all t > ty, 
then VY is either closed or approaches a closed 
trajectory along a spiral with the passage of 
time. 

We illustrate this in Figure 88. Here D 
consists of two closed curves IT, and T, 
and the circular domain between them. 
With each boundary point of D we associate 
a vector 


Viz, y) = X (@, y)it+Y G@, y)j. 


Then, if a trajectory [ that emerges at the 
initial moment ¢ = t, from a boundary 
point, enters D, and remains there at all 
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moments ¢2t1,, then, according to the 
above criterion, it will along a_ spiral 
approach a closed trajectory I, that lies 
entirely in D. The curve I, must surround 
a singular point of the differential system, 
point P, not lying in D. 

The differential system (173) provides 
a simple example illustrating the appli- 
cation of the above criterion in finding limit 
cycles. Indeed, system (173) has only one 
singular point, O(0, 0), and, therefore, 
the domain D lying between the circles 
with radii r = 1/2 and r = 2 contains no 
singular points. The first equation in (177) 
implies that dr/dt is positive on the inner 
circle and negative on the outer. Vector V, 
which is associated with the points on the 
boundary of D, is always directed into D. 
This means that the circular domain lying 
between circles with radii r = 1/2 and 
r = 2 must contain a closed trajectory of the 
differential system (473). Such a _ closed 
trajectory does indeed exist, it is a circle of 
radius r = 1. 

Note, however, that great difficulties are 
generally encountered in a system of the 
(422) type when we wish to realize prac- 
tically the Poincaré-Bendixson criterion, 
since no general methods exist for building 
the appropriate domains and, therefore, 
success depends both on the type of system 
and on the experience of the researcher. At 
the same time we must bear in mind that 
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finding the conditions in which no limit 
cycles exist is no less important than estab- 
lishing criteria for their existence. In this 
respect the most widespread condition is the 
Dulac criterion: if there exists a function 
B (a, y) that is continuous together with its 
first partial derivatives and is such that in 
a simply connected domain D of the phase 
plane the sum 


0 (BX) 0d (BY) 
Ox Oy 


is a function of fixed sign,* then no limit 
cycles of the differential system (122) can 
exist in D. At B (a, y) =1 the criterion 
transforms into the Bendizson criterion. 

[f we turn to the differential equation 
(156), which constitutes a differential model 
for describing an adiabatic one-dimensional 
flow of a perfect gas of constant specific 
heat ratio through a nozzle with drag, then 
for this equation we have 


X (a, y)=(1—y) F (2), 
¥ (x, y)=4y (41+ = y) (ray—F' (a). 
If we put 

B(x, y)= {4yP (2) [4+ 2 = y |} 


* That is, positive definite or negative definile, 
16* 
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we find that 


0 ee ) Ae vq 
ar a “— F (2) > 


and, ice (156) has no closed integral 
curves. 

Let us discuss one more concept that can 
be employed to establish the existence of 
limit cycles. The concept is that of the 
index of u singular point. 

Let [T be a simple closed curve (i.e. a 
curve without self-intersections) that is 
not necessarily a phase trajectory of system 
(122), lies in the phase plane, and does not 
pass through the singular pvints of this 
system. Then, if P (zx, y) is a point of IT, 
the vector 


with i and j the unit vectors directed along 
the Cartesian axes, is a nonzero vector and, 
hence, is characterized by a certain direction 
specified by an angle @ (Figure 89). 
point P (x, y) moves along I, say, counter- 
clockwise and completes a full cycle, vector 
V performs an integral number of cycles 
in the process, that is, angle 0 acquires an 
increment AQ = 2nmn, with n a_ positive 
integer, a negative integer, or zero. This 
number n is said to be the indez of the closed 
curve I’ (or the index of cycle I). 

If we begin to contract T in such a manner 
that under this deformation IT does not 
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pass through singular points of the given 
vector field, the index of the cycle must, 
on the one hand, vary continuously and, 
on the other, remain an integer. This means 
that under continuous deformation of the 
curve the index of the cycle does not change. 
This property leads to the notion of the 
index of a singular point as the index of a 
simple closed curve surrounding the sin- 
gular point. 

The index has the following properties: 

(1) the index of a closed trajectory of the 
differential system (122) is equal to +1, 

(2) the index of a closed curve surrounding 
several singular points is equal to the sum 
of the indices of these points, and 

(3) the index of a closed curve encompassing 
only ordinary points is zero. 
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This implies, for one thing, that since 
the index of a closed trajectory of system 
(122) is always +14, a closed trajectory must 
enclose either one singular point with an 
index +1 or several singular points with 
a net index equal to +1. This fact is often 
used to prove the absence of limit cycles. 

The index of a singular point is calculated 
by the formula 


n=1+ 5S", (179) 


where e is the number of elliptic sectors, 
and h is the number of hyperbolic sectors. 
For practical purposes the following simple 
method may be suggested. Suppose that LZ 
is a cycle that does not pass through sin- 
gular points of (422) and is such that any 
trajectory of (122) has no more than a finite 
number of points common to L. The tra- 
jectories may intersect ZL or touch it. In 
the latter case only exterior points of con- 
tact (type A) or interior points of contact 
(type B) are taken into account, while the 
C-type points, points of inflection, are not 
(see Figure 90). We can still use formula 
(179) to calculate the index of a singular 
point, but now e is the number of interior 
points of contact and h the number of 
exterior points of contact of trajectories 
of (122) with cycle ZL. 

Figure 91 depicts singular points with 


= 
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indices 0, +2, +3, +1, and —2, re- 
spectively. 

Earlier it was noted that constructing 
a complete picture of the behavior of the 
phase trajectories of the differential system 
(122) is facilitated by introducing the point 
at infinity via transformations (159). To- 
pological considerations provide a very 
general theorem which states that when 
a continuous vector field with a finite num- 
ber of singular points is specified on a 
sphere, the net index of the points is +2. 
Thus, if the net index of all the singular 
points of a differential system (possessing 
a finite number of such points) that lie in 
a finite phase-plane domain is distinct 
from +2, the point at infinity must be a 
singular point with a nonzero index. 

But if instead of inversion we employ 
homogeneous coordinates, the net index of 
all the singular points is already +41. That 
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this is so can be seen from the fact that if 
the plane is projected onto a sphere with 
the center of projection placed at the center 
of the sphere, two points on the sphere cor- 
respond to a single point on the projective 
plane, and the circumference of the great 
circle parallel to the plane corresponds to 
a straight line at infinity. 

If we turn to Eq. (168), which describes 
adiabatic one-dimensional flow of a perfect 
gas through a rotating tube of uniform cross 
section, then, as shown in Figure 85, for the 
singular point at infinity we have e = 2 
and h = 0. It follows from this that the 
index of this point is +2 and does not 
depend on the values of the constants in 
Eq. (168). This implies, for one thing, that 
the net index of finite singular points is 
zero. It can also be shown that, depending 
on whether the straight line specified by 
the equation my — Gz? = 0 has two points, 
one double point, or not a single point of 
intersection with the parabola fixed by the 
equation 1— py + qG*x? = 0 (which, in turn, 
is equivalent to G being greater than Go, 
equal to Gy, or Jess than Go, with Go = 
2mq'/?/p), the following combinations of 
finite singular points occur. 

(a) G>>G,. A saddle point and a nodal 
point, 

(b) G = G,. A higher-order singular point 
with two hyperbolic (hk = 2) sectors and 
two parabolic (e = 0). 
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(c) G<G,. Singular points are absent. 
We see that in three cases the net index 
is zero, as it should be. 


2.144 Periodic Modes in Electric 
Circuits 


We will show how limit cycles emerge in 
a dynatron oscillator (Figure 92). An anal- 
ysis of the operation of such an oscillator 
leads to what is known as the Van der Pol 
equation. Although phenomena linked with 
generation of limit cycles can be illustrated 
by examples from mechanics, biology, and 
economics,’ we will show how such phe- 
nomena emerge in the study of electric 
circuits. 

Figure 92 depicts schematically a dyna- 
tron oscillator, with the i,-v, characteristic 
shown by a solid curve in Figure 93. Here 
i, is the current and v, the voltage in the 


Fig. 92 
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screen-grid tube, the dynatron. The circuit 
consists of resistance r, inductance JZ, 
and capacitance C connected in parallel 
and known as the tank circuit, which is 
connected in series with the screen-grid tube. 
The real circuit can be replaced in this case 
with an equivalent circuit shown schematic- 
ally in Figure 94. The characteristic of the 
tube may be approximated with a third- 
degree polynomial i = av + yv*, which is 
shown by a dashed curve in Figure 93. 
Here i and v stand for the coordinates in 
a system whose origin is shifted to the 
point of inflection O. As follows from 
Figure 93, 


a>0, y>O0. 

In accordance with one of Kirchhoft’s laws, 
i+i, +i, +ie =0, 

with i, =v/r, L (di,/dt) =v, and ig = 


Cy. As a result of simple manipulations we 
arrive at the following differential equation: 


3 
v+(F +34+43 oe ve) V+ T= 0. 
If we now put 


1 


| 3 
otaean Gab eae 


we can write the previous equation in the 
form 


v + (a + bv) v + wo =0. 
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Fig. 94 


This differential equation is known as the 
Van der Pol equation. If we introduce the 
transformations 


ease Ney 


we can associate with the Van der Pol 
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equation the following first-order differen- 
tial equation: 


dy ss (a+ bv*) y+ojv _ Y(v, y) 
‘dv y is V (v, y) oe) 


The only finite singular point of this 
equation is the origin, and 


J* (0, 0) = wo? > 0, 
D (v, y) = —(a + bv*). 


Since 6 >0O, we assume that a>O and 
conclude that divergence D does not re- 
verse its sign and, hence, Eq. (180) can 
have no closed integral curves. We, there- 
fore, will consider only the case where a 
is negative, that is, a < —1/r. This im- 
plies that D (0, 0) = —a > 0 and, hence, 
the singular point is either a nodal point 
or a focal point. If we now consider the 
differential system corresponding to the 
differential equation (180), that is, 


dv 


dt = Y(v, y), 2 tee CL y); 


we see that as ¢ grows, the representative 
point moves along a trajectory toward the 
singular point. Thus, a trajectory that 
emerges from the point at infinity cannot 
reach the singular point at the origin no 
matter what the value of ¢ including the 
case where ¢ = + oo. This implies that if 
we can prove that a trajectory originating 
at the point at infinity resembles a spiral 
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winding onto the origin as t > oo, this 
will guarantee that at least one limit cycle 
exists. 

The existence proof for such a cycle can 
be obtained either numerically or analyti- 
cally. A numerical method suggested by the 
Dutch physicist and mathematician Van 
der Pol (1889-1959) consists in building 
a trajectory that originates at a point 
positioned at a great distance from the 
origin and in checking whether this tra- 
jectory possesses the above-mentioned prop- 
erty. Such a procedure yields an approxi- 
mation to the limit cycle, but it can be 
employed only in the case where concrete 
numerical values are known. 

Below we give a proof procedure* based 
on analytical considerations and the in- 
vestigation of the properties of singular 
points at infinity. Here, in contrast to the 
method by which Eq. (168) was studied, 
we employ the more convenient (in the 
present problem) honugeneous-cvordinate 
transformations 


y= E/2, ens: (181) 


The straight line at infinity is fixed by the 
equation z = 0. To reduce the number of 


* See the paper by J. Kestin and S.K. Zaremba, 
“Geometrical methods in the analysis of ordinary 
differential equations,” Appl. Sci. Res. B3: 144 
189 (1953). 
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variables to two, we first exclude the point 


§ = z= 0. This can be done if we assume 
that € = 1. Then 


De 42. yf yi. 
In this case the differential system asso- 
ciated with Eq. (480) is transformed to 


dz dy (az? +-b)y + w222 
a ee ee ee 


« 


It is convenient to introduce a new para- 
meter, 0, in the following manner: 


di = 22 dd. (182) 


Then this system of equations can be writ- 
ten as 


ce — (az? + b) n—wf22*— 22? = Y, 

183) 
d ( 
w= = 


Note, first, that the straight line at infinity, 
z =U, constitutes a trajectory of the dif- 
ferential system (183) and as ¢ grows (hence, 
as 8 grows) the representative point moves 
along this trajectory toward the only sin- 
gular point z = yn = 0. The characteristic 
equation, which in this case has the form 
—byz = 0, specifies the regular exceptional 
direction z = 0, with two repulsive regions 
corresponding to this direction, which can 
be established. by studying the appropriate 
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Fig. 95 


vector field (Figure 95). The second excep- 
tional direction yn =O is singular and, 
therefore, additional reasoning is required. 

The locus of points ¥ = 0 is a curve that 
touches the straight line » =O at the 
origin and passes through the second and 
third quadrants (Figure 95), which makes 
it possible to fix three directions, J, JJ, 
and J/I, on different sides of the symmetry 
axis z = Q. The region lying between the 
curve fixed by the equation ¥Y =O and 
the axis ny = 0 is topologically equivalent 
to two repulsive regions. Thus, on each side 
of the straight line y = O there is at least 
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one trajectory that touches the straight 
line at the origin. 

If we now consider the differential equa- 
— 


at — |-— weit t=! (0, z), 


we note that for small values of | 7 | 


Of w? 
5 — (1-4) <0 
in the second quadrant between curve Y = 
0 and axis y = 0. Hence, if we take two 
phase trajectories with the same value of 1, 
it is easy to see that the representative 
points moving along these trajectories will 
also move apart as z decreases. This means 
that on each side of the straight line y = 0 
there can be only one phase trajectory that 
touches this straight line at the origin 
and belongs to the region considered. It is 
also clear that the qualitative picture is 
symmetric about the axis z = 0. More than 
that, since for small values of z and positive 
nm we have 


[W/Z | S w?/| yz |, 


there can be no curves that touch the axis 
7 = 0 at the origin and pass through the 
first or fourth quadrant. Such curves are 
also absent from the region that lies to the 
left of the curve fixed by the equation 
& =O since in this case ¥ is positive and 


17—0770 
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2 has the same sign as z (Figure 95, 
arrows IY). This reasoning suggests that 
the singular point y» = z= 0 is a saddle 
point. The two phase trajectories that reach 
this point as t + -+-oo are segments of the 
straight line at infinity, z = O, that link the 
previous point with the point § = z= 0. 
The other two phase trajectorics reach the 
saddle point as t --— oo. 

In our reasoning we did not discuss the 
point § = z = 0. To conclude our investi- 
gation, let us put 7 = 1 in (181). Then the 
differential system associated with the dif- 
ferential equation (180) can be written as 
follows: 


SEP +E (az? + DE Off24) = 


Sy = 2 (a2? + OE + otc?) = Q, 


a 


where the variable 6 is defined in (182). The 
point € = z = 0 proves to be singular, and 
the characteristic equation here is 2? = 0. 
Hence, the exceptional direction z = 0 is 
singular, too. The curve fixed by the equa- 
tion P (€, z) = 0, (Figure 96) touches the 
axis 2 = O at the origin and has a cuspidal 
point there, while the curve fixed by the 
equation Q (€, z) = 0 has a multiple point 
at the origin. By studying the signs of P 
and Q we can establish the pattern of the 
vector field (Figure 96) and the phase-tra- 
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Fig. 97 


jectory pattern in the neighborhood of the 
singular point € = z = 0 (Figure 97). 

We note, for one thing, that far from the 
axis &°= 0 there are no phase trajectories 
that enter the singular point from the right. 


17* 
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This follows from the fact that in the first 
and fourth quadrants 


| 2/E | > | Q/P |. 


The above reasoning shows that the Van 
der Pol equation has no phase trajectories 
that tend to infinity as f grows but has an 
infinitude of phase curves that leave infinity 
as ¢ grows. This proves the existence of at 
least one limit cycle for the Van der Pol 
differential equation. 


2.15 Curves Without Contact 


In fairly simple cases the complete pattern 
of the integral curves of a given differential 
equation or, which is the same, the phase- 
trajectory pattern of the corresponding 
differential system is determined by the 
type of singular points and closed integral 
curves (phase trajectories), if the latter 
exist. Sometimes the qualitative picture 
can be constructed if, in addition to estab- 
lishing the types of singular points, we can 
find the curves, separatrices, that link the 
Singular points. Unfortunately, no general 
methods exist for doing this. Hence, it is 
expedient in qualitative integration to 
employ the so-called curves without con- 
tact. The reader will recall that a curve or 
an arc of a curve with a continuous tangent 
is said to be a curve (arc) without contact 
if it touches the vector (X, Y) specified 
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by the differential system (122) nowhere. 
The definition implies that vector (X, Y) 
must point in one direction on the entire 
curve. Thus, a curve without contact may 
be intersected by the phase trajectories of 
(122) only in one direction when ¢ grows 
and in the opposite direction when ¢ di- 
minishes. Therefore, knowing the respective 
curve may provide information about the 
pattern of a particular part of a phase tra- 
jectory. 

Various auxiliary inequalities can be 
used in qualitative integration of differen- 
tial equations. For example, if two differ- 
ential equations are known, say, 


dy dy 
Ge i y), de 8 (, y) 


and it is known that f(z, y)< g (z, y) 
in a domain D, then, denoting by y, (z) 
a solution of the first equation such that 
Y1 (%o) = Yo. With (%, Yo) ED, and by 
y, (x) a solution of the second equation 
with the same initial data, we can prove 
that y, (7) < yy, (z) for x >, in D. But 
if in D we have the strict inequality 
f (x, y)<g (@, y), then y; (x) < yp (2) for 
x>z, in D and the curve y = y, (2) is 
a curve without contact. 

As an example let us consider the differ- 
ential equation (168). We have already 
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shown that if G > G,, the equation allows 
for two finite singular points, a saddle 
point and a nodal point, which are points 
of intersection of the straight line my — 
G’x =0 and the parabola 41— py + 
qgG*x? = 0. The segments of the straight line 
and parabola that link these two finite sin- 
gular points are curves without contact and 
they specify a region of the plane which we 
denote hy A. If we put 


X (x, y) = 1 — py + 9G2’, 
Y (x, y) = 2y (my — Ga), 


it is easy to see that vector (X, Y) is di- 
rected on the boundary of A_ outward 
except at the singular points. Hence, if the 
representative point emerges from any inte- 
rior point of A and moves along an integral 
curve with ¢ decreasing, it cannot leave A 
without passing through one of the singular 
points. But since inside A we _ have 
X (x, y) <0, the singular point that is 
attractive is the nodal point. 

Finding the slopes of the exceptional 
directions for the saddle point, we can see 
that one of the exceptional straight lines 
passes through A. This implies that an 
integral curve that is tangent to this 
straight line at the saddle point must enter 
A and then proceed to the nodal point. 
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2.146 The Topographical System 
of Curves. The Contact Curve 


Earlier we noted that in the qualitative 
solution of differential equations it is expe- 
dient to use curves without contact. It must 
be noted, however, that there is no general 
method of building such curves that would 
be applicable in every case. Hence, the 
importance of various particular methods 
and approaches in solving this problem. 
One approach is linked to the selection of 
an appropriate topographical system of 
curves. Here a topographical system of curves 
defined by an equation ® (x, y) = C, with 
C a real-valued parameter, is understood to 
be a family of nonintersecting, enveloping 
each other, and continuously differentiable 
simple closed curves that completely fill a 
doubly connected domain G in the phase 
plane.* 

If the topographical system is selected 
in such a manner that each value of para- 
meter C is associated with a unique curve, 
then, assuming for the sake of definiteness 
that the curve corresponding to a definite 
value of C envelopes all the curves with 
smaller values of C (which means that the 
“size” of the curves grows with C) and that 
no singular points of the corresponding 


* Other definitions of a topographical system also 
exist in the mathematical literature, but essentially 
they differ little from the one given here. 
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differential system lie on the curves be- 
longing to the topographical system, we 
arrive at the following conclusion. If we 
consider the function ® (q@ (£), w (t)) where 
x = gp (t) and y = > (¢) are the parametric 
equations of the trajectories of the differen- 
tial system of the (122) type and calculate 
the derivative with respect to ¢, that is, 


d@ t), OO (2, 
wee tp (2) =_ Cs y) X (x, y) 
0 

{ew YY (x, y), 

then the curves (cycles in our case) without 
contact are the curves of the topographical 
system on which the derivative d@/dt 
is a function of fixed sign. If, in addition, 
d@/dt is positive on a certain curve be- 
longing to the topographical system, then 
such a curve, being a cycle without contact, 
possesses the property that all phase tra- 
jectories intersecting it leave, with the 
passage of time ¢, the finite region bounded 
by it. And if d@/dt is negative, all the 
phase trajectories enter this region. This 
implies, for one thing, that if in an annular 
region completely filled with curves belong- 
ing to the topographical system the deriv- 
ative d@/dt is a function of fixed sign, 
the region cannot contain any closed tra- 
jectories, say, the limit cycles of the dif- 
ferential system. Limit cycles may exist 
only in the annular regions where the 
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derivative d@/dt is a function with alter- 
nating signs. 

To illustrate this reasoning, we consider 
the differential system 


d d 
sr=y—eta, Sba—r—yty?, (184) 


which has only one singular point in the 
finite part of the plane, O (0, 0), a stable 
focal point. For the topographical system 
of curves we select the family of concentric 
circles centered at point O (0, 0), that is, 
the family of curves specified by the equa- 
tion 2? -+ y? =C, with C a positive pa- 
rameter. In our case the derivative d@/dt 
is given by the relationship 


< = —2(x2+ y?) 4+ 2 (xt4+-y4), 


which in polar coordinates x = r cos 0 and 
y=rsin@ assumes the form 


c= — 2r? + 2r* (cos! 8 + sin‘ 0) = — 2r? 


+ 2r* ( S++ cos 40 ) : 


Noting that the greatest and smallest values 
of the expression in parentheses are 1 and 
1/2, respectively, we arrive a the following 
conclusion. For r greater than ) 2 the value 
of the derivative d@/dt is positive while 
for r less than unity it is negative. On the 
basis of the Poincaré-Bendixson criterion 


18-0770 
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we conclude (if we replace ¢ with —# in 
system (184)) that the annular region bound- 
ed by the circles 2? + y27=1 and 2? + 
2 == 2 contains a limit cycle of system 
(184), and this limit cycle is unique. 

To prove the uniqueness it is sufficient 
to employ the Dulac criterion for a doubly 
connected domain: if there exists a function 
B (x, y) that is continuous together with 
its first partial derivatives and is such that 
in adoubly connected domain G belonging 
to the domain of definition of system (122) 
the function 


a (BX) (BX) 4 9 (BY) ey) 


is of fixed sign, then in G there can be no 
more than one simple closed curve consisting 
of trajectories of (122) and containing within 
it the interior boundary of G. 

In our case, selecting B (zx, y) =1 for 
system (184), we find that 0X/dr + 
OY /dy = 3 (a? + y?) — 2 and, as can easily 
be seen, in the annular region with bound- 
aries 2? -+ y2=1 and 2?+y? =2 the 
expression 3 (x? + y’?) — 2 always retains 
its sign. Allowing now for the type of the 
singular point O (0, 0), we conclude that 
the cycle is an unstable limit cycle. 

To return to the problem of building 
curves without contact in the general case, 
let us examine a somewhat different way 
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of using the topographical system of 
curves. This approach is based on the 
notion of a contact curve. Introduction of 
this concept can be explained by the fact 
that if the derivative d@/dé vanishes on 
a set of points of the phase plane, this set 
constitutes the locus of points at which the 
trajectories of the differential system touch 
the curves belonging to the topographical 
system. Indeed, the slope of the tangent to 
a trajectory of the differential system is 
Y/X and the slope of the tangent to a curve 
belongirg to the topographical system is 
—(00 dzx)/(AD/dy). Thus, when 


av aD y, | 

ie Atay Y= 0; (185) 
the slopes coincide, or 

Y aD | a® 


Xo el oy 


Hence, the lucus of points at which the 
trajectories of the differential system of the 
(122) type touch the curves belonging to the 
topographical system (defined by the equa- 
tion ® (2, y) = C) is called a contact curve. 
The equation of a contact curve has the 
(185) form. Of particular interest here is 
the case where the topographical system 
can be selected in such a manner that either 
the contact curve itself or a real branch of 
this curve proves to be a simple closed 
curve. Then the topographical system con- 


1) 
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Fig. 98 


tains the “greatest” and “smallest” curves 
that are touched either by a contact curve 
or by a real branch of this curve (in Fig- 
ure 98 the contact curve is depicted by a 
dashed curve). If the derivative d@/dt 
is nonpositive (nonnegative) on the “great- 
est” curve specifed by the equation 
® (x, y) = C, and nonnegative (nonposi- 
tive) on the “smallest” curve specified by 
the equation ® (z, y) = C,, then the annu- 
lar region bounded by these curves contains 
at least one limit cycle of the differential 
system studied. 

Specifically, in the last example (see 
(184)), the contact curve specified by the 
equation 27+ y7=a2'-+ y* is a closed 
curve (Figure 99), which enables us to find 
the “greatest” and “smallest” curves in the 
selected topographical system that are 
touched by the contact curve. Indeed, if 


= 


Ch. 2. Qualitative Methods 269 


Fig. 99 


we note that in polar coordinates the equa- 
tion of the contact curve has the form 7? = 
(cost @ + sin* 6)-!, we can easily find the 
parametric representation of this curve: 


om cos 9 y= sin 8 
V cos? 0-F sin‘ 6 ’ Y cos! 6+ sin 0 


Allowing now for the fact that the greatest 
and smallest values of r® are, respectively, 
2 and 1, we conclude that the “greatest” 
and “smallest” curves of the topographical 
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system 2? + y? =C are the circles x? -- 
yy>=1 and 2?+ y*=2,_ respectively. 
These circles, as we already know, specify an 
annular region containing the limit cycle 
of system (484). 


2.17 The Divergence of a Vector 
Field and Limit Cycles 


Returning once more to the general case, 
we note that in building the topographical 
system of curves we can sometimes employ 
the right-hand side of system (122). It 
has been found that the differential system 
may be such that the equation 


aX oY 
where A is a real parameter, specifies the 
topographical system of curves. 

For instance, if we turn to the differential 
system (184), Eq. (186) assumes the form 
3 (x? + y?) — 2 =A. Then, assuming that 
A = (A+ 2)/3, with 4 € (—2, +00), we 
arrive at the topographical system of curves 
xz* 4+- y* = A used above. 

A remark is in order. If the “greatest” and 
“smallest” curves merge, that is, the contact 
curve or a real branch of this curve coin- 
cides with one of the curves belonging to 
the topographical system, such a curve is 
a trajectory of the differential system, 


Ch. 2. Qualitative Methods 274 


Let us consider, for example, the differ- 
ential system (473): 


d 
= yt (i-w—y, 


ae x+y (A— x2 — y?). 


For this system 

Ox 

je tay tA) 

and Eq. (186) assumes the form 2 — 
4 (x? + y”) =i. If instead of parameter A 
we introduce a new parameter A = 
(2 — A)/4, with % € (—oo, 2), then the 
latter equation, which assumes the form 
x? + y* = A, defines the topographical 
system of curves. The contact curve in this 
case is given by the equation (z? + y?) x 
(x? -+ y* — 1) = 0 and, as we see, its real 
branch x? + y* = 1 coincides with one of 
the curves belonging to the topographical 
system. This branch, as shown on pp. 237- 
239, proves to be a limit cycle of system 
(173). 

The last two examples, of course, illus- 
trate, a particular situation. At the same 
time the very idea of building a topograph- 
ical system of curves by employing the 
concept of divergence proves to be fruitful 
and leads to results of a general nature. 
We will not dwell any further on these re- 
sults, but only note that the following 
fact may serve as justification for what 
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has been said: if a differential system of the 
(122) type has a limit cycle L, there exists a 
real constant X and a positive and continu- 
ously differentiable in the phase plane func- 
tion B (x, y) such that the equation 


a Be) 4. 3 (BY) ey 


A, 


specifies a curve that has a finite real branch 
coinciding with cycle L. 

Let us consider, for example, the differ- 
ential system 


oe 22 — 223 -+ x*y — 3xy?-- y%, 


oy — a+ xy + xy’, 


which has a limit cycle specified by the 
equation x? + y? = 1. For this system the 
divergence of the vector field is 


a a 


We can easily verify that there is not a 
single real A for which the equation 2 — 
ox? — 3y* = Xd specifies a trajectory of the 
initial differential system. But if we take 
a function B (a, y) = 32? — 4azy + Ty? + 
3, then 


8 (BX) “ 6 (BY) 
dy 


Oz ol cal 
x (— 2372+ 16zy — 25y? — 20) — 14, 
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and the equation 
a(BX) , 0(BY) _ 


Ox | Oy usG 


specifies a curve whose finite real branch 
zg? + y» —1=0 proves to be the limit 
cycle. 

In conclusion we note that in studying 
concrete differential models it often proves 
expedient to employ methods not discussed 
in this book. Everything depends on the 
complexity of the differential model, on 
how deep the appropriate mathematical 
tools are developed, and, of course, on the 
erudition and experience of the researcher. 
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Appendix 1. 
Derivatives of Elementary 


Functions 
Function Derivative 
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Function Derivative 
4 
tan x a= sec? x 
COS? x 
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sin*z 
sin x 
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sinh x cosh zx 
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Function 


cosh7! x 
tanh7z 


coth—! x 


Appendix 2. 
Basic Integrals 


Power-Law Functions 


\ o"de= = (n = 1) 


Trigonometric Functions 


sinzdz= —coszx 
cosxdzr=sinz 


tan x dz= — In | cosz | 


! 
\ cot zdz=In | sing | 
| oe 
| ars 


Fractional Rational Functions 


{ita toe (2) 
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\ a =— coth=! ( =| 
1 r— 


=, ne (|x| >a) 
Exponential Functions 
{ iowa 

\ 0 dae 


Hyperbolic Functions 


neds =eoshe 
\ cosh zdz=sinhz 
tanh zdz= In | cosh | 


coth zdz=In | sinhz | 


d 
(= cok a tanh z 
\ =a = —cothz 


<ahts 


Irrational Functions 


| Taam tieh (=} 
ear Es: (=) =n on en 
\ Ta a cosh (= | =In | (2+-V z?—a?)| 


To the Reader 


Mir Publishers would be grateful for 
your comments on the content, transla- 
tion and design of this book. We would 
also be pleased to receive any other 
suggestions you may wish to make. 
Our address is: 

Mir Publishers 

2 Pervy Rizhsky Pereulok 

I-110, GSP, Moscow, 129820 

USSR 


